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I Abstract. In this paper we present a compilation of results from our most advanced neutron star merger simu- 

^ . lations. Special aspects of these models were refered to in earlier publications (Ruffert & Janka 1999, Janka et 

' al. 1999), but a description of the employed numerical procedures and a more complete overview over a large 

^SJ , number of computed models are given here. The three-dimensional hydrodynamic simulations were done with a 

^Nj ■ code based on the Piecewise Parabolic Method (PPM), which solves the discretized conservation laws for mass, 

' momentum, energy and, in addition, for the electron lepton number in an Eulerian frame of reference. Up to 

, five levels of nested cartesian grids ensure higher numerical resolution (about 0.6 km) around the center of mass 

• while the evolution is followed in a large computational volume (side length between 300 and 400 km). The sim- 

' ulations are basically Newtonian, but gravitational-wave emission and the corresponding back-reaction on the 

, hydrodynamic flow are taken into account. The use of a physical nuclear equation of state allows us to follow 

^ l' the thermodynamic history of the stellar medium and to compute the energy and lepton number loss due to 

the emission of neutrinos. The computed models differ concerning the neutron star masses and mass ratios, the 
neutron star spins, the numerical resolution expressed by the cell size of the finest grid and the number of grid 
^ levels, and the calculation of the temperature from the solution of the entropy equation instead of the energy 

, equation. The models were evaluated for the corresponding gravitational-wave and neutrino emission and the 

• • ■ mass loss which occurs during the dynamical phase of the merging. The results can serve for comparison with 

, ' smoothed particle hydrodynamics (SPH) simulations. In addition, they define a reference point for future models 

, with a better treatment of general relativity and with improvements of the complex input physics. Our simulations 

}_( ' show that the details of the gravitational-wave emission are still sensitive to the numerical resolution, even in 

our highest-quality calculations. The amount of mass which can be ejected from neutron star mergers depends 
strongly on the angular momentum of the system. Our results do not support the initial conditions of temperature 
and proton-to-nucleon ratio needed according to recent work for producing a solar r-process pattern for nuclei 
around and above the A « 130 peak. The improved models confirm our previous conclusion that gamma-ray 
bursts are not powered by neutrino emission during the dynamical phase of the merging of two neutron stars. 

Key words, stars: neutron — binaries: close — hydrodynamics — gravitational waves — nuclear reactions, nucle- 
osynthesis, abundances — elementary particles: neutrinos 



1. Introduction 

The binary pulsar PSR 1913-^16 (Hulse & Taylor 1975) is 
the most famous example for a binary system containing 
two neutron stars, among another ~ 1000 of such systems 
expected to exist in our Galaxy. High-precision measure- 
ments show that the change in time of the orbital parame- 
ters of PSR 1913-1-16 is consistent with expectations from 
the theory of general relativity, which predicts the emis- 
sion of gravitational waves and a continuous decrease of 
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the orbital separation. Therefore, these systems have a fi- 
nite lifetime of typically hundreds of millions up to billions 
of years. As the two stars spiral in towards each other, the 
evolution accelerates because the gravitational- wave emis- 
sion rises strongly with decreasing distance. When the or- 
bital separation has shrunk to only a few stellar radii, 
the system has become a strong source of gravitational 
waves with a frequency around 100 Hz. It will end its life 
within milliseconds in the final, catastrophic merging of 
the two neutron stars, emitting a powerful outburst of 
gravitational radiation which carries important informa- 
tion about the properties of the merging stars, the dynam- 
ics of the coalescence, and the remnant left behind. 



2 



M. Ruffert and H.-Th. Janka: Coalescing neutron stars - a step towards physical models 



With an estimated rate of about 10 ^ events per 
year per galaxy (e.g., sec the recent numbers in Buhk 
et al. 1999, Fryer et al. 1999, Kalogera & Lorimer 2000, 
and references therein) neutron star mergers are among 
the most frequent and most promising candidates for 
gravitational-wave emission which is strong enough to 
be measurable by the upcoming interferometric experi- 
ments in the U.S. (LIGO), Europe (GEO600, VIRGO), 
and Japan (TAMA) (Thorne 1995). Theoretical models 
and wave templates, however, are needed to help filter 
out the weak signals from disturbing background noise. 
Gravitational waves from neutron star mergers could be 
one of the most fruitful ways to learn about the internal 
properties of neutron stars. 

Merging neutron stars are also considered as possible 
sources of at least the subclass of short and hard cos- 
mic gamma-ray bursts, especially if the merger remnant 
collapses to a black hole on a dynamical timescale (for re- 
cent discussions and model calculations, see, e.g., Popham 
et al. 1999, Ruffert & Janka 1999). Coincident detections 
of gravitational waves and gamma rays would be a con- 
vincing observational confirmation of this hypothesis and 
might in fact be the only possibility to identify the cen- 
tral engine of a gamma-ray burst unequivocally. The X- 
ray satellite HETE-2, which was launched in Fall 2000, is 
hoped to bring a similar breakthrough in the observation 
of short bursts as the BeppoSAX satellite did in case of 
the long ones. 

The energy of the relativistically expanding fireball 
or jet, which finally produces the observable gamma-ray 
burst, can be provided by the annihilation of neutrino- 
antineutrino pairs (Paczynski 1991, Mcszaros & Rees 
1992, Woosley 1993a) or possibly by magnetohydrody- 
namical processes (Blandford & Znajek 1977, Meszaros 
and Rees 1997). In the former case, the gravitational bind- 
ing energy of accreted disk matter is tapped, in the latter 
case the rotational energy of the central black hole could 
be converted into kinetic energy of the outflow. If neutrino 
processes are supposed to power the gamma-ray burst phe- 
nomenon, very high neutrino luminosities are needed, of 
magnitude similar as those from core-collapse supernovae. 
The rate of neutron star mergers, however, is much smaller 
(by a factor of 100-10000) than the Galactic supernova 
rate. This practically excludes them as detectable sources 
of thermal neutrinos in the MeV energy range, because the 
signals are too faint to be measurable from extragalactic 
distances. Dissipative processes in the rclativistic outflow, 
which are considered to produce the gamma-ray burst, 
may also lead to the generation of high-energy or even ul- 
tra high-energy neutrinos (Paczynski & Xu 1994, Waxman 
& Bahcall 1997, 2000). Such neutrinos might be seen in 
future km^-scale experiments like ICECUBE, which is cur- 
rently under construction in the Antarctica. However, they 
do not carry much specific information about the origin of 
the relativistically moving particles and it is therefore not 
very likely that they can yield much evidence about the 
nature of the central engine that powers the gamma-ray 
burst. 



Neutron-rich matter, which is ejected from the system 
during the dynamical phase of the merging, was suggested 
as a possible site for the rapid neutron capture process (r- 
process) to produce heavy nuclei beyond the iron group 
(Lattimer et al. 1974, 1976; Hilf et al. 1974; Eichler et 
al. 1989; Meyer 1989). This problem has gained new in- 
terest recently (Rosswog et al. 1999, 2000; Preiburghaus 
et al. 2000). The possible contribution to the Galactic r- 
process material is estimated from the gas mass that gets 
unbound during the violent last stages of the coalescence. 
The nuclear reactions in decompressed neutron star mat- 
ter depend sensitively on the initial conditions (neutron 
excess, composition, temperature, density), the dynam- 
ical and, in particular, thermal history of the material, 
and the influence of beta-decays and corresponding neu- 
trino losses. All of these issues are so far not well under 
control in theoretical models, and therefore hydrodynamic 
simulations of neutron star mergers have not (yet?) been 
able to yield conclusive results. 

These questions have been the motivation for a large 
number of investigations of the spiral-in phase and the ul- 
timate merging of neutron stars. Analytic studies and el- 
lipsoidal treatments concentrated on the effects of viscous 
dissipation for the heating and the rotation of the stars 
(Kochanek 1992, Bildsten & Cutler 1992, Lai 1994), the 
final instability of the mass transfer near the tidal radius 
(e.g., Bildsten & Cutler 1992; Lai et al. 1994a,b; Taniguchi 
& Nakamura 1996; Lai & Wiseman 1996; Lombardi et 
al. 1997; Baumgarte 2001) and the deformed equilib- 
rium structure and tidal lag of the binary configuration 
prior to the dynamical interaction (Lai & Shapiro 1995). 
Hydrodynamical simulations of the coalescence were per- 
formed for Newtonian gravity with SPH codes (e.g., Rasio 
& Shapiro 1992, 1994, 1995; Centrella & McMillan 1993; 
Zhuge et al. 1994, 1996; Davies et al. 1994; Rosswog et 
al. 1999, 2000) and with grid-based methods (e.g., Oohara 
& Nakamura 1990, Nakamura & Oohara 1991), partly in- 
cluding special treatments of the gravitational-wave emis- 
sion and their back-reaction on the flow by adding the 
corresponding post-Newtonian terms to the equations of 
hydrodynamics (e.g., Ruffert et al. 1996, 1997a; Ruffert et 
al. 1997b). More recently progress has been achieved in a 
wider use of the post-Newtonian approximation (Shibata 
et al. 1998, Ayal et al. 2001, Faber & Rasio 2000, Faber 
et al. 2001) and considerable advances were made towards 
general relativistic treatments (Oohara & Nakamura 1999, 
Shibata 1999; Shibata & Uryii 2000, 2001). 

A spectacular result was obtained by Mathews & 
Wilson (1997, and references therein) who found that rel- 
ativistic effects lead to a compression of the two nciitron 
stars during the late stages of the spiral-in and therefore 
to their gravitational collapse to black holes prior to the 
merging. This effect contradicts Newtonian models where 
tidal stretching reduces the density of the stars as they 
get closer. Analytic considerations confirm the Newtonian 
behavior also for the post-Newtonian case (Thorne 1998; 
Baumgarte et al. 1998a, b), and more recent simulations 
by the Wilson group (Marronetti et al. 1999) as well as 
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general relativistic hydrodynamic models by other groups 
(Bonazzola et al. 1999, Shibata et al. 1998) were not able 
to reproduce the result of Mathews & Wilson (1997). The 
latter was recognized to be due to an error in the approx- 
imation scheme to full general relativity (Flanagan 1999). 
In any case, pre-merging collapse of the neutron stars is 
a speculative option only if the nuclear equation of state 
is extraordinarily soft and the neutron stars are already 
very close to the maximum mass for stable single neutron 
stars. 

The majority of the simulations by other groups was 
done with simple microphysics, in particular with a poly- 
tropic law P = Kp^ for the equation of state (EoS) of the 
neutron star matter. This is a fair approach when one is 
mainly interested in the calculation of the gravitational- 
wave emission, which is associated with the motion of the 
bulk of the mass. It offers the advantage that the influence 
of the stiffness of the EoS, which determines the mass- 
radius relation of the neutron stars and the amount of 
compression which occurs during the final plunge, can be 
easily studied by choosing different values for the adiabatic 
index V. 

Several years ago we started to compute merger models 
with a more elaborate treatment of the EoS of the neutron 
star matter, using the physical description by Lattimer & 
Swesty (1991), which enabled us to follow the thermody- 
namics of the gas and to include a treatment of the neu- 
trino production and emission from the heated neutron 
stars (Ruffert et al. 1996, 1997a; Ruffert & Janka 1998, 
1999; Janka et al. 1999). Our main aims were the investi- 
gation of the relevance for gamma-ray burst scenarios, in 
particular for those where the neutrino emission had been 
suggested to provide the energy for the relativistic gamma- 
ray burst fireball via neutrino-antineutrino annihilation. 
Also the amount of mass ejection during the dynamical 
interaction and the properties of the ejected matter de- 
pend on the EoS, which cannot be descibed by one simple 
polytropic law in both the low-density and high-density 
regimes. 

After publication of our first papers (Ruffert et 
al. 1996, 1997a), we changed our code considerably and, 
in particular, we improved many features which had in- 
fluence on the results of our simulations. For example, 
we introduced nested grids to get a higher resolution of 
the neutron stars and at the same time to use a larger 
computational volume. In addition, we extended the EoS 
table to higher temperatures and lower densities. The lat- 
ter allowed us to reduce the density of the dilute medium 
that has to be assumed around the neutron stars on the 
Eulerian grid. Since the heat capacity of cold, degenerate 
matter is very small, minor numerical noise in the internal 
energy had induced larger errors in the temperature. We 
therefore also implemented an entropy equation, because 
the entropy is numerically less problematic for calculat- 
ing the temperature. Besides these improvements, we also 
covered a wider range of scenarios, e.g., added models with 
opposite directions of the neutron star spins and with dif- 
ferent neutron star masses as well as different mass ratios. 



All of our later publications referred to data of models 
which were computed with the improved version of the 
code. So did the simulations of the black hole accretion 
in Ruffert & Janka (1999) start from an initial model of 
the new generation of calculations, and also in the tables 
of Janka et al. (1999) data of new neutron star merger 
models were listed. So far, however, we published only very 
specific aspects of these new models and did not present 
our results in detail. This is the purpose of the present 
publication. 

In Sect. H we will give a technical description of the 
code changes and improvements, in Sect. ^ a list of com- 
puted models, in Sect. |4|we shall present the main results 
for the new models, and in Sect. H we shall discuss the 
implications and draw conclusions. 

2. Changes and improvements of the numerics 

The three-dimensional Eulerian hydrodynamical conser- 
vation laws for mass, momentum and energy for non- 
viscous flow are integrated explicitly in time on a carte- 
sian grid, using a finite-volume scheme which is based 
on the Piecewise ParaboHc Method (PPM) of Colella & 
Woodward (1984). The code is basically Newtonian, but 
includes the post-Newtonian terms that account for the 
local effects of gravitational-wave emission (the volume 
integral of these local terms reproduces the quadrupole 
approximation) and the corresponding back-reaction on 
the hydrodynamical flow according to the formulation by 
Blanchet et al. (1990) (for details, see Ruffert et al. 1996). 
The Poisson equations for the gravitational potential and 
two additional potentials used in the back-reaction terms 
are solved by fast Fourier transforms. For the described 
neutron star merger simulations, the thermodynamics of 
the stellar medium are described by a tabular version of 
the Lattimer & Swesty (1991) EoS, which allows us to take 
into account the source terms for energy and lepton num- 
ber loss due to neutrino production. These source terms 
are calculated with a neutrino trapping scheme which eval- 
uates the production rates of neutrinos and antineutrinos 
of all flavors. The trapping scheme takes into account the 
optical depth at any location inside the star to reduce 
the neutrino release when the diffusion timescale becomes 
long. Since the emission of electron neutrinos and antineu- 
trinos can change the electron lepton number of the stellar 
medium, we also solve a continuity equation for this quan- 
tity. A detailed discussion of the equations and a more 
complete explanation of the numerical methods can be 
found in Ruffert et al. (1996, 1997a). 

The numerical code and the input physics described 
in Ruffert et al. (1996, 1997a) have been changed and im- 
proved in a number of aspects. 

(i) Grid: 

The cartesian grid for the simulations presented there had 
a side length of only 82 km. Therefore a signiflcant amount 
of mass was swept off the grid although it would not have 
become unbound. For this reason we enlarged the compu- 
tational volume to 328-656 km by introducing up to flve 
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Fig. 1. Illustration 
of four levels of 
nested grids in two 
(instead of three) 
spatial dimensions, 
each level with 64 
cells in every di- 
rection, covering a 
computational vol- 
ume with side length 
of 328 km. The inner- 
most two grid levels 
are enlarged and the 
initial positions of 
the two neutron stars 
in the orbital plane 
are indicated 



levels of nested grids (Fig. |^). This ensures good numeri- 
cal resolution near the grid center so that the two neutron 
stars are represented in the new calculations with a sim- 
ilar accuracy as the best resolved Model A128 in Ruffert 
et al. (1996, 1997a). The larger computational volume al- 
lows us to trace the matter with high angular momentum, 
which is flung out to large distances but returns and is 
added to the cloud of gas that surrounds the compact 
merger remnant. 

(ii) EoS table: 

Another problem in the previous calculations was that 
under certain extreme conditions (e.g., extreme heating as 
in case of Model C64) the boundaries of the EoS table were 
hit. Therefore we expanded the temperature range in a 
new EoS table to cover values from 10~^ MeV to 100 MeV, 
and reduced the minimum density to 5 x lO^gcm"^. 

(in) Environmental density: 

Since the simulations are done in an Eulerian reference 
frame, the density outside the neutron stars cannot be set 
to zero but only to a small value, i.e., a value small com- 
pared to the average density in the stellar interior. In order 
not to influence the dynamical behavior of the ejected stel- 
lar fragments by gas that is swept up in the surroundings, 
the density of this ambient medium should be chosen as 
low as possible. Therefore we used a value of 10* g cm~^ 
in the new calculations instead of lO^gcm"'^ previously. 
The total mass of this dilute gas on the grid is therefore 
significantly less than 1O~^M0. Since only an extremely 
small fraction of this gas interacts with the neutron star 
matter, the dynamical effect of this medium is negligible. 

(iv) Temperature determination and entropy equation: 
The employed hydrodynamics code solves the energy 
equation for the total specific energy, which is the sum 
of the specific internal and kinetic energies. Compared to 
integrating an equation for the internal energy alone, this 
has the advantage that without gravity the energy equa- 



tion is written in a fully conservative way. The effects of 
gravity are included by source terms. It has, however, the 
disadvantage that the internal energy, from which the tem- 
perature is determined, has to be calculated as the differ- 
ence of total and kinetic energies. If both the latter ener- 
gies are large, a small value of the internal energy can be 
significantly affected by numerical noise, and the tempera- 
ture determination can become inaccurate. This is usually 
not a problem, but does become problematic in the case 
of extremely degenerate neutron star matter, which has 
a very small heat capacity. Only minor variations of the 
internal energy can then lead to much larger changes of 
the temperature. For these reasons the temperature is a 
quantity which is very sensitive to numerical deficiencies. 
As a consequence, it was not possible to start the simu- 
lations in Ruffert et al. (1996, 1997a) with cold neutron 
stars. A stable temperature evolution could be obtained 
when the thermal energy density in the initial configu- 
ration was assumed to be two per cent of the degeneracy 
energy density. With this prescription, the central temper- 
ature of the neutron stars was several MeV and the surface 
temperature about half an MeV initially. Nevertheless, the 
corresponding neutrino luminosities were negligibly small, 
and the total thermal energy was so tiny compared to the 
other energies (internal, gravitational and kinetic) that its 
effect on the dynamics was unimportant. 

In order to avoid these problems in the temperature 
determination and to have the freedom of starting with 
colder neutron stars, we decided to use the entropy in- 
stead of the energy density for temperature calculations. 
The entropy evolution was followed by a separate entropy 
equation which was integrated in addition to the hydro- 
dynamic conservation laws. In an Eulerian frame of ref- 
erence, the corresponding continuity equation is (making 
use of Einstein's summation convention) 

d d(^srL\^v^'j 

— {s nb) H ^— = S'y + 5*811 + S^is ■ (1) 
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Here s means the matter entropy per nucleon, rib — p/u 
the baryon number density {u being the atomic mass unit, 
p the rest mass density) , and the cartesian components 
of the velocity vector. The divergence of the entropy flux 
was written in its cartesian form (summation is applied 
when an index appears twice), and the source terms on 
the right hand side are the rates of change of the entropy 
density due to neutrino production, shock dissipation, and 
shear and bulk viscosity effects, respectively. 

The entropy generation rate per unit volume by neu- 
trino processes is (e.g., Cooperstein 1988) 



1983, Landau & Lifschitz 1991) 



Se 



- S-L {ipe +'4'p- V'n) 



(2) 



where Se and S'l are the effective energy loss rate and 
the effective lepton number source term as defined in 
Appendix B of Ruffert et al. (1996), A:bT is the temper- 
ature in MeV and the tp^s denote the degeneracy param- 
eters (chemical potentials divided by the temperature) of 
electrons, protons and neutrons (including the rest mass 
energies) . The neutrino source terms as given by Ruffert et 
al. (1996) are only evaluated for temperatures above about 
0.5 MeV. Below this threshold they are switched off, be- 
cause the assumptions employed in their calculation are 
not valid any more. 

The entropy generation rate per unit volume by shocks 
is given according to the tensor formalism of Tscharnutcr 
& Winkler (1979) as 



•S'sh — 



Ql4 



(3) 



with the mixed tensor Q™ of the viscous pressure given 

by 



Pp 



3 dx^ J 



dx^ 



< 



(4) 



otherwise 



and the mixed tensor e™ of the symmetrized gradient of 
the velocity field defined by 



dvi 



dxi dx^ 



(5) 



where Einstein's summation convention is used (dv^ /dx'^ 
therefore means the divergence of the velocity vector in 
cartesian coordinates), and Sf^ is the mixed unity ten- 
sor. The characteristic length / is of the order of the local 
width of the grid, I = f-Ax. We calibrated the proportion- 
ality factor / such that the entropy jump across a shock 
as calculated with our hydrodynamics code is reproduced 
by the entropy generation according to Eqs. (|1)-(|). We 
found best agreement for the choice of / = 1.8. 

The entropy generation rate per unit volume due to 
shear and bulk viscosity can be written (using again 
the summation convention) as (e.g., Shapiro & Teukolsky 



5'vis — 



dv£\ 

dx^ J 



dx^ 



+ c 



dx'^ 



(6) 



where 77 is the dynamic shear viscosity coefficient and C 
the bulk viscosity coefficient. Solving the Euler equations 
for an ideal fluid, shear viscosity effects are only caused by 
the numerical viscosity of our hydrodynamics code, which 
we describe by the ansatz 77 = apvAx. With a typical grid 
resolution Ax between 10** cm and 10^ cm one empirically 
finds values for a between 5 x 10~^ and 5 x 10"'^ (Janka et 
al. 1999). We used a representative number of a = 2x 10^^ 
in the simulations discussed below. 

In hot neutron star matter with the proton frac- 
tion exceeding a critical lower limit, bulk viscosity can 
be strongly enhanced by the direct URCA processes of 
electron neutrino and antineutrino production and ab- 
sorption (Haensel & Schaeffer 1992). For matter com- 
posed of neutrons (n), protons {p) and electrons (e) with 
trapped neutrinos (but with no trapped lepton-number 
excess, i.e., if ipe + ''Pp ^ ^ 1) one can write an ap- 
proximate expression for the bulk viscosity coefficient as 
C 102'*(ypp/po)i/3 g/(cms) with Yp = Ye being the pro- 
ton fraction and po w 2.5 x 10^"' g/cm^ the density of 
normal nuclear matter (Haensel & Schaeffer 1992, Sawyer 
1980, van den Horn & van Weert 1981). Although the nu- 
merical factor was taken somewhat larger than estimated 
by Haensel & Schaeffer (1992), we found the entropy gen- 
eration rate associated with the bulk viscosity term to be 
negligibly small. 

The entropy loss rate due to neutrino emission, ex- 
pressed by the source term S'j^, is tiny initially, but be- 
comes (globally, i.e. as an integral over all grid cells) com- 
parable to the (positive) shear viscosity term ^vis after 
several milliseconds, when the neutron stars have merged 
to a hot, rapidly and differentially spinning object. Earlier 
than this, in particular prior to the merging, the time is 
too short for shear viscosity to raise the entropy, and the 
temperature is too low for neutrinos to make any effect. 
When the neutron stars begin to touch (roughly half a 
millisecond after the simulations were started), the shock 
dissipation term Sgh becomes clearly the dominant one 
globally, about twice to twenty times bigger than S'vis- 

Using Eq. (|l|) for evolving the entropy, an updated 
value of the temperature is obtained in a predictor-cor- 
rector step which ensures second order accuracy for the 
time integration. In a first step one evaluates the entropy 
source terms with the old temperature, then solves Eq. (|l]) 
to get an estimate for the new entropy and thus for the 
new temperature, and then solves Eq. (^ a second time 
with source terms computed with an average value of the 
old and estimated new temperature. 

The temperature thus obtained from the entropy equa- 
tion is used to calculate the neutrino source terms in the 
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hydrodynamic conservation laws of mass, momentum, en- 
ergy, and Icpton number, which still describe the evolu- 
tion of the stellar fluid. This is not fully consistent and 
should therefore not be considered as the necessarily bet- 
ter treatment. Instead, it is meant as an alternative ap- 
proach which allows one to test the uncertainties associ- 
ated with the temperature determination and the corre- 
sponding effects of neutrinos. 
(v) Neutrino treM,tm,ent: 

Let us conclude this section with a few remarks about the 
treatment of neutrinos. The energy radiated in neutrinos 
during the computed evolution (about 10 ms) is typically 
more than an order of magnitude smaller than the total 
energy emitted in gravitational waves. Whereas neutrino 
emission is very small during the first five milliseconds, it 
increases later and dominates the energy loss in the second 
half of the computed evolution. For the dynamical phase 
of the merging process (which lasts only a few millisec- 
onds after the start of the simulations), neutrino source 
terms are therefore insignificant. Of course, the treatment 
of neutrino effects by using a trapping scheme (Ruffert et 
al. 1996, 1997a) is a strong simplification of the true prob- 
lem. An exact treatment would require three-dimensional, 
time-dependent (general relativistic) transport of neutri- 
nos and antineutrinos in a moving neutron star medium, 
which is extremely optically thick in dense and hot re- 
gions and transparent near the stellar surface and in re- 
gions where the matter is cold with temperatures below 
about 1 MeV. Solving this problem is currently not feasi- 
ble, but we are convinced that our trapping scheme is a 
good first approach which yields an approximate descrip- 
tion of effects connected with the production of neutrinos, 
and a reasonably good estimate for the total luminosity 
of neutrinos. 

Neutrinos stream off freely as soon as they are pro- 
duced in transparent matter, thus changing the quantity 
(i.e., energy, entropy or lepton number) of the stellar 
medium according to the rate given by the neutrino source 
term Sv'- 



dt 



(7) 



In contrast, the dominant mode of energy, entropy and 
lepton number loss is by neutrino diffusion when neutrinos 
are in equilibrium with the stellar medium at high optical 
depths. This is expressed by the equation 



dt 



(8) 



The description used in the trapping scheme is guided by 
these cases. It calculates the loss terms for energy (and 
entropy) and lepton number from the local neutrino emis- 
sion rates in the optically thin regime, and from the rate 
of diffusive depletion of the equilibrium densities of neu- 
trino number and energy in the optically thick regime. 
The transition between both limits is done with a smooth 
interpolation based on the local diffusion timescale, which 
is estimated from the optical depth to the stellar surface 



(see Appendix B in Ruffert et al. 1996). In the equilibrium 

regime, howevcir. we neglect the contributions of neutrinos 
to the internal energy density and to the entropy density, 
or, in other words, we neglect that neutrinos in equilib- 
rium contribute to the heat capacity of the stellar matter 
(an effect which would affect the computed temperature). 
A similar approximation is also applied for the (electron) 
lepton number because we solve a continuity equation only 
for the net electron number density (number density of 
electrons minus number density of positrons) and neglect 
the lepton numbers carried by electron neutrinos and an- 
tineutrinos. These approximations can be justified by the 
fact that the chemical potential of the electrons is much 
larger than that of neutrinos. Therefore the net lepton 
number of electron neutrinos minus antineutrinos is very 
small, and the summed energy densities of neutrinos and 
antineutrinos of all flavors are always less than ~ 10% of 
the energy density of electrons (plus positrons, if present) 
plus photons, and even smaller compared to the total in- 
ternal energy density of the gas, which also includes the 
contributions of nucleons and nuclei. For the same reason, 
we also made no effort to take into account the neutrino 
pressure (or momentum transfer to the stellar medium 
when neutrinos start to decouple). 

The use of a trapping scheme means that neutrino 
diffusion or propagation relative to the stellar medium 
and neutrino advection along with the fluid motion arc 
ignored. Therefore the corresponding energy and lepton 
number (and entropy) transport inside the star are disre- 
garded. These effects are likely not to be very important 
within the few milliseconds of the computed evolution, be- 
cause the diffusion timescale in the dense, hot interior of 
neutron stars is several seconds and neutrino lepton num- 
ber and energy are small compared to the medium quanti- 
ties. Even in the less dense region between neutrinosphere 
and nuclear core the changes of energy and lepton num- 
ber due to local processes (which are evaluated with the 
trapping scheme) should dominate the transfer of energy 
and Icptons by neutrinos which stream from one region of 
the star to another. Therefore we think that the trapping 
scheme is suitable to account for the main effects which 
play a role on millisecond timescales. If the simulations 
were continued for a longer time, tens or hundreds of mil- 
liseconds, the transport effects would, of course, have to 
be included. 

The trapping scheme fails to yield a suitable approxi- 
mation when and where neutrino heating is stronger than 
neutrino cooling. This is the case outside of the neutrino- 
sphere. Here high-energy neutrinos transfer energy to the 
cooler stellar gas via scattering and absorption reactions. 
This energy deposition causes an outflow of baryonic mass, 
the so-called neutrino-driven wind, which has been inves- 
tigated in some detail for cooling proto-neutron stars in 
type-II supernovae (Duncan et al. 1986, Woosley 1993b, 
Qian & Woosley 1996). It is a major disadvantage that 
the trapping scheme does not allow one to study this in- 
teresting phenomenon in the context of merging neutron 
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Table 1. Computed models with their characterizing parameters and results. The models differ concerning the chosen 

initial spins of the neutron stars. In several cases the temperature evolution was followed by using an entropy equation, 
starting from "high" ("ent wa") or low initial temperature ("ent co"). Also a special choice of the grid is indicated 
by the model name and in the column with remarks. N is the number of grid zones per dimension in the orbital 
plane, g the number of levels of the nested grid, L the size of the largest grid, I the size of the smallest zone. Mi and 
M2 are the masses of the two neutron stars, ao is their initial center-to-center distance, /cbTo the initial maximum 
temperature (reached at the center of the neutron stars), fsim the computed period of time of the evolution, Mp<ii 
is the gas mass with a density below 10^^ g/cm'^ at the end of the simulation, Md the gas mass with specific angular 
momentum larger than the Keplerian angular momentum at a radius equal to three Schwarzschild radii of the merger 
remnant, Mg the mass of the gas which leaves the grid, Mu the gas mass which can become unbound, and ksTrnx 
the maximum temperature reached on the grid during the simulation. The > signs indicate that the numbers are still 
changing when the simulations were stopped. 



model remark spin N g L I Mi M2 ao fee To tsim Mp<ii Ma Mg Mu ksTm: 

km km M© Mq km MeV ms Mg/lOO Mg/lOO Mg/lOO Mq/100 MeV 



A32 




none 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


5.44 


10. 


4.5 


6.1 


0.78 


0.03 


53. 


A64 




none 


64 


4 


328 


0.64 


1.6 


1.6 


42. 


5.47 


10. 


5.8 


8.0 


1.78 


0.23 


39. 


B32 




solid 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


6.76 


10. 


6.9 


21. 


8.3 


1.6 


39. 


B32w 


ent wa 


solid 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


5.44 


10. 


6.5 


19. 


9.6 


2.2 


30. 


B32w' 


ent'wa 


solid 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


5.44 


10. 


6.1 




9.7 


2.2 


49. 


B32c 


ent CO 


solid 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


0.05 


10. 


6.5 


19. 


9.6 


2.1 


29. 


B32^c 


ent CO 


solid 


32 


5 


656 


1.28 


1.6 


1.6 


42. 


0.05 


24. 


14. 


18. 


4.5 


1.9 


30. 


B64 




solid 


64 


4 


328 


0.64 


1.6 


1.6 


42. 


6.25 


10. 


5.7 


25. 


9.2 


2.4 


39. 


B64c 


ent CO 


solid 


64 


4 


328 


0.64 


1.6 


1.6 


42. 


0.05 


10. 


5.5 


25. 


9.1 


2.4 


40. 


C32c 


ent CO 


anti 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


0.05 


10. 


4.2 


1.7 


0.13 


0.005 


58. 


C64c 


ent CO 


anti 


64 


4 


328 


0.64 


1.6 


1.6 


42. 


0.05 


10. 


5.7 


6.0 


0.35 


0.0085 


69. 


C128c 


ent CO 


anti 


128 


4 


328 


0.32 


1.6 


1.6 


42. 


0.05 


2. 


>0.3 


>1. 


>0. 


>0. 


78. 


032 




oppo 


32 


4 


328 


1.28 


1.6 


1.6 


42. 


0.05 


10. 


4.0 


9.0 


1.57 


0.36 


66. 


064 




oppo 


64 


4 


328 


0.64 


1.6 


1.6 


42. 


0.05 


10. 


3.4 


8.4 


1.19 


0.19 


89. 


S32 




solid 


32 


4 


328 


1.28 


1.2 


1.2 


42. 


4.68 


11.6 


6.0 


16. 


7.9 


2.4 


32. 


S64 




solid 


64 


4 


328 


0.64 


1.2 


1.2 


42. 


4.71 


10. 


6.5 


23. 


7.5 


2.0 


35. 


D32 




solid 


32 


4 


400 


1.56 


1.8 


1.2 


46. 


7.16 


10. 


7.1 


14. 


8.8 


2.7 


39. 


D64 




solid 


64 


4 


400 


0.78 


1.8 


1.2 


46. 


7.19 


13. 


7.0 


13. 


9.3 


3.8 


35. 


D64^ 


3 grids 


solid 


64 


3 


400 


1.56 


1.8 


1.2 


46. 


7.16 


5. 


>2.7 


22. 


>0.8 


>0.4 


>33. 


D128^ 


3 grids 


solid 


128 


3 


400 


0.78 


1.8 


1.2 


46. 


7.19 


2. 


>0.2 


>5. 


>0. 


>0. 


>12. 



stars. We shall try to remove this deficiency in future im- 
provements of our code. 

3. New simulations 

In this section we shall describe the different models and 
will review the most important results of our simulations. 

Binaries with different baryonic masses of the neutron 
stars, 1.2, 1.6, and 1.8 M0, were considered, with differ- 
ent mass ratios (1:1 and 1:1.5), and with different neutron 
star spins, where we distinguish the four cases of initially 
irrotational systems ("none"), synchronous (or 'tidally 
locked') rotation ("soHd"), counterrotation ("anti"), and 
opposite spin directions ( "oppo" ) of the two neutron stars. 
The angular velocity due to the spins, which is added to 
the orbital motion, is equal to the angular velocity of the 
orbit at the chosen initial distance of the two stars, but 
the spin directions are varied between the four cases (sec 
also Ruffert et al. 1996). In addition, we computed mod- 
els with different grids, i.e., with different zone sizes of 



the finest grid (and thus different numerical resolution) as 
well as different numbers of grid levels. Finally, we used 
the entropy equation to follow the temperature history of 
our models, and started the computations with different 
temperatures of the neutron stars, "warm" models with 
a central temperature of 5 7 MeV, satisfying the require- 
ment that the thermal energy at the beginning is two per 
cent of the internal energy, and "cold" models with a con- 
stant initial temperature of only 0.05 MeV. 

3.1. Initial conditions 

The initial configuration consists of two spherically sym- 
metric neutrons stars, constructed as hydrostatic equi- 
librium solutions for Newtonian gravity and the EoS of 
Lattimer & Swesty (1991). The temperature of cold and 
warm models is chosen as mentioned above, and the profile 
of the electron fraction, Y(.{r), corresponds to the equi- 
librium state of cold, deleptonized, neutrino-transparent 
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neutron stars (i.e., the chemical potential of electron neu- 
trinos vanishes throughout the star). The 1.6 Mq neutron 
star has a radius of about 15 km, the 1.2 Mq star is slightly 
smaller and the 1.8 Mq star a little bigger. This scaling of 
mass and radius indicates that the effective adiabatic in- 
dex of the stellar matter is larger than two (see Shapiro & 
Teukolsky 1983). The initial center-to-center distance oo 
of the two stars was set to a value such that the stars could 
finish about one full revolution before the final plunge oc- 
curred at the radius of tidal instability. Circular orbits 
were assumed for all models. The initial orbital veloc- 
ity was chosen according to the inspiraling motion of two 
point masses Mi and M2 at separation ap in response to 
their emission of gravitational waves. From the quadrupole 
formula the angular velocity and the radial velocity can 
be calculated as (Cutler & Flanagan 1994): 



lG{Mi+M2) 



64 1 
5 Oq 



Vr = a = ~ — ^^ MiM2(Mi + M2) 



(9) 
(10) 



More details about the treatment of the initial conditions 
can be fomid in Ruffert et al. (1996). 

3.2. Models 

Table |l| contains a list of computed models with their 
specific characteristics. Basically we distinguish models of 
types A,B,C,0,S and D. Models A,B,C and O are our 
standard cases with two equal neutron stars with baryonic 
masses of 1.6 Mq and initially irrotational (spin: "none"), 
corotational ("solid"), counterrotational ("anti") and op- 
posite directions ("oppo"), respectively, for the spins in 
the initial state. The S-models are computed with smaller, 
1.2 Mq neutron stars, and the D-models with two different 
neutron stars of 1.2 Mg and 1.8 Mg. (Note that Model V64 
in Janka et al. 1999 is identical to Model C64c of this 
work.) 

The names of the models carry information also about 
the number of zones of each grid level (32, 64, or 128), 
and about the number of grid levels, indicated by a super- 
script, if it is not the standard value of 4. Since we take 
the extension of the grid in the z-direction perpendicular 
to the orbital plane to be only half as big as in the orbital 
plane and, in addition, assume equatorial symmetry, the 
grids have 32x32x8 or 64x64x16 or 128x128x32 zones, 
respectively. A higher grid level has only half the reso- 
lution (twice the zone size) of the level below. With an 
equatorial length and width of the computational volume 
between 328 and 656 km, the smallest zones have a side 
length between 0.32 km and 1.56 km. For fixed computa- 
tional volume a smaller number of grid levels implies that 
the individual grids are bigger, and avoids that the two 
neutron stars cross the boundary of the innermost grid 
during spiral in, as they do for our standard case of four 
levels of nested grids (see Fig. |l|) . Reducing the number of 



grid levels therefore allows one to test the corresponding 
numerical effects. 

Use of the entropy instead of the energy to calculate 
the temperature (which is then used for computing the 
neutrino emission) is indicated by an extension of the 
model name. The letter "w" marks the "warm" models, 
"c" the "cold" ones. Since our stars have to be embed- 
ded by a medium with finite density on the Eulerian grid, 
the motion of the stars through this medium creates a 
shock wave at the stellar surfaces, where the tempera- 
ture increases in a narrow ring of grid zones. This effect 
is energetically umimportant, and also does not lead to a 
significant increase of the neutrino luminosity as long as 
the volume and the density of the heated medium stay 
low. However it is unphysical and can lead to an inflation 
of the surface layers of the neutron stars. Therefore we 
try to reduce it by localizing grid cells with temperature 
spikes in the shocked region below a certain density and 
resetting the temperature there to an average value of the 
surrounding grid zones. The influence of this manipulation 
has to be tested. On the one hand we did this by comput- 
ing models with increasingly better resolution (where the 
disturbing effects occured in a smaller volume of space), 
on the other hand we changed the procedure for detecting 
the zones with unphysical temperatures. A corresponding 
test model is B32w'. 

4. Results 

In this section the results of our new simulations will be 
described as far as important differences to our previ- 
ously published models showed up, or interesting effects 
occurred in dependence of the parameters varied between 
the models. 

4.1. Dynamical evolution and gravitational-wave 
emission 

Figures |^-|^ show the hydrodynamical results for five of 
the models listed in Table |l|. The displayed cases in- 
clude Model A64 with two equal, initially nonrotating 
neutron stars (Fig. ||), Model B64 with two equal, ini- 
tially corotating stars (Fig. ||), Model C64 with two equal, 
initially counterrotating stars (Fig. ^), Model 064 with 
equal neutron stars with the spins pointing in opposite 
directions (Fig. and Model D64 with neutron stars 
of different masses in initially locked rotation (Fig. 
Figure ^ provides cuts perpendicular to the orbital planes 
for Models B64, 064, and D64, whereas the other plots 
display the density and temperature distributions in the 
orbital plane. 

Compared to our previously published simulations 
(Ruffert et al. 1996, 1997a), Models A64, B64, and C64 
were computed with a finer resolution around the grid 
center and simultaneously a significantly larger volume 
by employing nested grids. Corresponding models with 32 
zones on the finest grid level have a similar resolution as 
the standard models in our preceding papers. In addition. 
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Model C64 was started with a lower initial temperature 
and the temperature history was followed by using the 
entropy equation. 

The displayed results reveal a significant dependence 
of the dynamical evolution on the neutron star masses 
and spins, fn the corotating case (B64) pronounced spiral 
arms form, which are much less strongly developed in the 
A and C models and are present only for the corotating 
component in Model 064. In case of the different neutron 
star masses (D64) the less massive star is stretched to a 
long, banana shaped object before it is finally wrapped 
around the more massive component. Prior to this, a siz- 
able amount of matter is flung out from the outward point- 
ing end of the smaller star. This material forms a spiral 
arm which expands away from the merger site. 

The initial total angular momentum of the binary sys- 
tems varies with the spins of the two components. Among 
the models with the same masses of the two neutron stars, 
the total angular momentum is largest for the B cases 
(neutron stars corotating with the orbit), smaller for the 
A (no neutron star spins) and O models (spins in opposite 
directions), and smallest for the C models (both neutron 
stars in counterrotation with the orbital motion) (Fig. 
During the computed evolution, about 10% of the angular 
momentum are removed by the emission of gravitational 
waves (Fig. ||). This leads to correlations between max- 
ima of the gravitational- wave luminosity and periods of a 
decrease of the angular momentum on the grid (compare 
Figs. 1^ and ^l]). More angular momentum is carried away 
by the matter that leaves the boundaries of the computa- 
tional grid. While the gravitational-wave emission peaks 
within the first 2 milliseconds, the effect due to the mass 
loss becomes important only later than 4-6 ms after the 
start of the simulations. In addition, the code does not 
conserve the angular momentum exactly. The violation 
depends mainly on the resolution of the finest grid level, 
where the bulk of the matter is located. For the mod- 
els with a cell size of 0.64 km on the central grid, less 
than 10% of the initial angular momentum were destroyed 
within about 10 ms of evolution. This, unfortunately, does 
not reach the excellent quality of the conservation of en- 
ergy by our code (Fig. ^ . 

The total angular momentum of the coalescing binary 
system determines the structure and the dynamical state 
of the very compact central body of the merger remnant, 
which contains the bulk (Jj 90%) of the system mass. 
The wobbling and ringing of this body after the merg- 
ing is sensitive to the fluid motions and thus to the value 
and distribution of angular momentum in its interior. The 
plots in Fig. |l^ show that the two stars in Model B64 fuse 
smoothly and form a centrally condensed body immedi- 
ately after the plunge, whereas in Model A64, and even 
stronger in Model C64, two density maxima can still be 
distinguished later. In particular in the counterrotating 
Model C64 the distance of these density maxima varies 
with time, indicating that the remnant remains in a per- 
turbed internal state. Correspondingly, the gravitational 
waveform of Model C64 (Fig. |l^ reveals a low-freqency 
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Fig. 8. Total angular momentum (component perpendic- 
ular to the orbital plane) of the matter on the grid (upper 
curves) and cumulative angular momentum loss by gravi- 
tational wave emission (lower curves) as functions of time 
for different models 

modulation on top of the high-frequency structure which 
is caused by the oscillation period of the compact remnant. 
This feature is essentially absent in Model B64. 

Before the final plunge, which is marked by the max- 
ima of the gravitational waveforms, the gravitational wave 
emission follows the chirp signal of two point masses rather 
accurately. We have not started our simulations with a 
configuration in rotational equilibrium and correspond- 
ingly deformed neutron stars. Therefore the neutron stars 
oscillate around their new quasi-equilibrium state before 
the merging, and the density maximum on the grid shows 
periodic modulations. Although these internal pulsations 
of the two stars are not damped by numerical viscosity 
(which expresses the non-dissipative quality of the code; 
Fig |l^) , they do not have any visible consequences for the 
gravitational waveforms. 

The deviation from the chirp signal becomes large 
when the orbital instability sets in and the neutron stars 
get strongly deformed and finally begin to touch each 
other. This moment contains very important information 
about the properties of the neutron stars, in particular 
about their mass-radius relation, which depends on the 
stiffness of the nuclear EoS. Similarly important and char- 
acteristic information about the binary system is carried 
by the post-merging signal, which is emitted if the mas- 
sive and compact merger remnant is stabilized by pressure 
or centrifugal forces and does not collapse to a black hole 
on a dynamical timescale. The collection of gravitational- 
wave amplitudes displayed in Fig. |l^ supports this argu- 
ment. The waveforms for the binary systems vary strongly 
with the spins of the neutron stars (compare Model O 
with Models A, B, and C) and with the mass ratio of 
the two stars (Model D). Of course, following the post- 
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Fig. 9. Total internal, kinetic, and gravitational potential energies of the matter on the grid as functions of time for 
Models B64 and C64c. The total energy is calculated by adding up these energies plus the energy that was emitted in 
gravitational waves. It should be conserved, because the gas that leaves the grid, has a total energy very close to zero 



Fig. 10. The separation of the density maxima and the values of the maximum density on the grid as functions of 
time for different models 



merging evolution definitely requires a general relativistic 

description of the problem, and quantitatively meaningful 
predictions of the gravitational-wave emission cannot be 
made on grounds of our basically Newtonian simulations. 
Relative differences between the models in dependence of 
the binary parameters should, however, also survive in a 
relativistic treatment. 



Comparison of the gravitational-wave amplitudes of 
Models A64, B64, and C64 with results displayed in 
Fig. 25 of Ruffert et al. (1996) reveals discrepancies for 
the post-merging signals. Because of the better numeri- 
cal resolution of our present simulations (the cell size on 
the finest grid is 0.64 km instead of 1.28 km for the old 
models), the numerical loss of angular momentum is sig- 
nificantly lower now. The influence of the resolution can 
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Fig. 11. Gravitational- wave luminosity and cumulative energy emitted in gravitational waves as functions of time for 
different models. The thin solid curves represent the emission of binary systems with two point masses instead of the 
neutron stars 



be directly tested by Models A32, B32, and C32, which 
have a factor of two larger cells on the innermost grid 
and thus have the same resolution as the old calculations. 
The structure of the wave amplitude for Models A32, B32, 
and C32 is therefore much more similar to our previously 
published results. 

The gravitational-wave luminosity (Fig. Ill]) peaks 
when the two neutron stars plunge into each other and the 
quadrupole moment of the binary changes most rapidly. 
For Models A, B, C, O this happens at roughly the same 
time, but the maximum luminosity of Models A and B 
is about twice as large as the one of the counterrotating 
Model C and of Model O where the neutron star spins are 
in opposite directions. Model S with smaller neutron stars 
and Model D with different neutron stars show signifi- 
cantly lower peak luminosities. For a short time, the lumi- 
nosity rises above the value of two orbiting point masses 
at the same distance. Because of the finite size of the neu- 
tron stars and the rapid development of a nearly spher- 
ical merger remnant, however, the luminosity decreases 
quickly after the maximum and never diverges as in case 
of the arbitrarily close approach of point masses. In the 
cases with strong, coherent fluid motion within the rem- 
nant, secondary and even tertiary maxima can occur. Such 
effects are absent in Model D and very weak in Model O. 

Our models are basically Newtonian, therefore the cal- 
culated gravitational waveforms are of limited usefulness 
as templates for measurements. Nevertheless, the simula- 
tions can serve for comparison with future, fully general 
relativistic models and will help one understanding the 
properties of their gravitational-wave signals. 

4.2. Dynamical mass loss 

Immediately after the merging of the neutron stars, tidal 
arms begin to reach out through the outer Lagrange points 



of the binary system. A little later, between 3 and 6 mil- 
liseconds after the start of the simulations, these spiral 
arms become inflated, because the matter expands with 
radial velocities up to one third of the speed of light. At 
that time the spiral arms reach the boundaries of the com- 
putational volume and mass flows off the grid in our sim- 
ulations. This phase coincides roughly with the moments 
shown in the middle panels of Figs. and can be rec- 
ognized in Fig. |l^ from the steep increase of the plotted 
curves. We monitor the fraction of the matter which ful- 
fills the condition that its total specific energy as the sum 
of gravitational, kinetic and internal energies is positive. 
This matter is considered to potentially become unbound. 

The amount of matter which is stripped off the tips of 
the spiral arms depends on the total angular momentum 
of the coalescing binary system. For solid-body type initial 
rotation we find the largest values: Roughly a tenth of a 
solar mass leaves the grid, and about 20-30% of this mat- 
ter might escape the system. The mass ejection is roughly 
a factor of ten smaller for the cases of irrotational neutron 
stars and neutron star spins in opposite directions, and 
another factor of ten smaller for the configuration with 
initially counterrotating stars (Table |l] and Fig. ^3|) . 

The gas mass Afg that is swept off the grid and the 
mass Mu that may get unbound exhibit only a rather 
weak dependence on the grid resolution, which can be 
seen by comparing models with different choices of grids 
m Table |l|. Of course, the size of the largest grid has a 
significant influence on Mg, but again does not affect 
very much (Model B32^c). 

Although the correlation of mass ejection with the an- 
gular momentum of the binary system appears plausible, 
and basically is consistent with the results of Rosswog et 
al. (2000), we cannot claim to have determined final val- 
ues for the dynamical mass loss. There are several factors 
of uncertainty which affect our numbers. The finite envi- 
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Fig. 13. The cumulative amount of gas that flows off the grid as function of time for different models (loft) and the 
amount of matter that is estimated to become unbound, because its total energy (as the sum of gravitational, kinetic, 
and internal energies) is positive 

Table 2. Gravitational-wave and neutrino emission properties^for all models. C, is the maximum gravitational-wave 
luminosity, E is the total energy emitted in gravitational waves, rh is the maximum amplitude of the gravitational waves 
scaled with the distance to the observer, r, L^^ is the electron neutrino luminosity after approaching a saturation level at 
about 8 ms, Lp^ is the corresponding electron antineutrino luminosity, and L^^ the luminosity of each individual species 
of heavy- lepton neutrinos, (ej/_.), (cp,) and (ei,^) are the mean energies of the emitted v^. z/g, and respectively. 
gives the total neutrino luminosity at the end of the simulation and E^jj denotes the integral rate of energy deposition 
by neutrino- antineutrino annihilation. 



model L Z rh L^^ Lo^ L^^ (c^) (e^^) (e,.^) E^a 

10=5 erg 10=2 -^qA 10=3 2££ 10==* ££S 10=3 £££ 10=3 2££ MeV MeV MeV 10=°2£S 
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A32 


2.3 


4.9 


8.6 


0.8 


2.0 


0.43 


4.5 


12. 


18. 


27. 




A64 


2.1 


5.2 


8.6 


0.9 


2.6 


0.37 


5.0 


11. 


17. 


27. 


90.5 


B32 


2.2 


3.8 


9.2 


0.6 


2.3 


0.45 


4.7 


12. 


17. 


25. 




B32w 


2.2 


3.5 


8.4 


0.65 


1.9 


0.23 


3.5 


13. 


15. 


25. 




B32w' 


2.2 


3.5 


8.4 


0.39 


1.3 


>0.08 


2.0 


11. 


15. 


22. 




B32c 


2.1 


3.5 


8.4 


>0.45 


1.5 


>0.09 


>2.3 


11. 


15. 


22. 




B32'5c 


2.1 


3.7 


8.3 


0.65 


1.0 


0.08 


1.9 


12. 


16. 


23. 




B64 


2.1 


3.7 


8.9 


0.6 


1.8 


0.22 


3.3 


12. 


17. 


25. 


70.0 


B64c 


2.1 


3.6 


8.5 


>0.62 


>1.6 


>0.07 


>2.5 


13. 


16. 


24. 




C32c 


0.85 


3.4 


6.0 


>1.0 


>2.2 


>0.2 


>4.0 


11. 


16. 


27. 




C64c 


1.2 


2.3 


6.0 


>1.1 


2.3 


0.13 


>4.0 


11. 


16. 


26. 




032 


1.3 


1.5 


6.9 


0.4 


0.9 


0.02 


1.4 


12. 


19. 


23. 




064 


1.2 


1.5 


6.9 


0.4 


1.1 


0.02 


1.6 


14. 


17. 


24. 




S32 


0.65 


1.4 


5.2 


0.4 


1.2 


0.18 


2.3 


11. 


16. 


25. 




S64 


0.70 


1.4 


5.5 


0.3 


0.9 


0.07 


1.5 


11. 


16. 


25. 




D32 


0.46 


0.96 


5.8 


0.4 


1.4 


0.19 


2.6 


12. 


17. 


24. 




D64 


0.37 


1.25 


5.5 


0.4 


1.0 


0.16 


2.0 


12. 


19. 


26. 





ronmental gas density, which we have to assume around 
the neutron stars for numerical reasons, has probably neg- 
ligible dynamical influence. More problematic is the de- 
creasing resolution of the nested grids towards the bound- 
aries of our computational volume. The largest uncer- 
tainty, however, is connected with the fact that we cannot 
follow the expanding matter until it moves ballistically. 



Therefore it is unclear how efficiently internal energy is fi- 
nally converted into bulk kinetic energy by hydrodynam- 
ical effects (i.e., PdV-work). Our results for the ejected 
mass can only be considered as best estimates, supported 
by the fact that the expansion of the tidal tails is mainly 
a kinematic effect and only to a minor degree driven by 
pressure forces. 
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Fig. 14. Electron fraction in the orbital plane at two different times for Model B64 (upper row) and Model D64. 
The contours are spaced linearly in steps of 0.02, beginning with a mininnun value of 0.02, which is also adopted for 
the very dilute environmental medium. The bold lines correspond to values of 0.02, 0.06, 0.1, 0.16, 0.2, and 0.3. The 
maximum values are 0.28 in the upper left plot, 0.18 in the upper right one, 0.35 in the lower left one, and 0.33 in the 
lower right one 



Matter which is ejected during the merging of binary 
neutron stars may have important implications for the 
production of heavy elements in our Galaxy. Dynamical 
events involving neutron stars (mergers, explosions) have 
been speculated to be possible sources of r-process nu- 
clei, both because of the high neutron densities present 
in the interior of neutron stars and because of the 
neutron-rich nuclei which exist in the crust below the 
surface of the cold neutron star (e.g., Hilf et al. 1974; 
Lattimer & Schramm 1974, 1976; Lattimer et al. 1977; 
SymbaUsty & Schramm 1982; Eichler et al. 1989; Meyer 
1989; Preiburghaus et al. 1999). The crust below a den- 
sity of about 2 X lO^^gcm"^ (see, e.g., Weber 1999) can 



contain up to several hundredths of a solar mass for a neu- 
tron star with a gravitational mass around 1.4 M0 (e.g., 
Akmal et al. 1998), although the exact value is somewhat 
EoS dependent. With an estimated merger rate of one per 
10^=*=^ years and a mean mass loss per event of ~ 10~^ 
solar masses, about 10-1000 solar masses of r-process ma- 
terial might have been produced by 10^ to 10^ such events 
over the history of our Galaxy. 

Quantitative predictions of the abundance yields are 
not easy to obtain. Ideally, the nuclear reactions and beta 
decays of neutron-rich isotopes should be followed along 
with the hydrodynamic flow of the expanding ejecta, be- 
cause heating by beta decays and neutrino losses can influ- 
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ence their thermal and chemical evolution. Moreover, it is 
essential to start the nucleosynthesis calculations with ap- 
propriate initial conditions. This concerns the initial tem- 
perature, initial electron fraction = ne/n\^, and the 
initial nuclear composition of the surface material of cold 
neutron stars. 

A first attempt for a nucleosynthetic evaluation of the 
merger ejecta was performed as a post-processing step of 
hydrodynamical data (Freiburghaus et al. 1999). Instead 
of being sclf-consistently determined as the consequence 
of neutrino reactions, the electron fraction was set to a 
chosen value. In addition, the initial temperature of the 
ejected neutron star matter was considered to be high 
enough (T ^ 5 x 10^ K) for the composition to be de- 
termined by nuclear statistical equilibrium (NSE). This is 
not necessarily a correct assumption, because the matter 
that is flung out in the tidal arms is not subject to any 
heating processes which could significantly increase the 
surface temperature of the originally cold neutron stars. 
In three-dimensional hydrodynamical simulations of neu- 
tron star mergers, however, the limited numerical resolu- 
tion does not allow one to accurately trace the thermal 
history of such cold matter, but numerical viscosity und 
numerical noise lead to artificial heating (for our code, 
such problems are discussed in Sects. |^ and ^). If the tem- 
perature stays low, which is most likely the case in the 
ejected and decompressed matter, the initial nuclear com- 
position of the neutron star crust is not erased by the 
onset of NSE. Instead, the neutron-rich nuclei start to 
undergo beta decays and the final distribution of stable 
or long-living nuclei is sensitive to the initial composition 
(e.g., Hilf et al. 1974, Lattimer et al. 1977, Meyer 1989, 
Sumiyoshi et al. 1998). 

In our simulations the electron fraction evolves as a 
result of the emission of electron neutrinos Ve and an- 
tineutrinos Vf.. Near the tips of the tidal tails, however, 
the neutrino producing reactions are not fast enough to 
cause a sizable change of the initial value of Yg on the 
short timescales of the dynamic expansion. This is so be- 
cause the temperature near the neutron star surface was 
chosen to be "only" around 1-2 MeV in the beginning (see 
Fig. 2 in Ruffert et al. 1996; such values of the tempera- 
ture are certainly unphysically high for cold initial stars, 
but they can be considered as low for the present discus- 
sion) and the temperature does not increase in the spiral 
arms during the subsequent evolution (Figs. ||-^. As a 
consequence, Yg can serve as a tracer quantity, and our re- 
sults indicate that only matter that was initially located 
close to the surfaces of the original stars is ejected dur- 
ing coalescence. This matter essentially retains its initial, 
very neutron-rich state with Ye significantly less than 0.1 
(which corresponds to the neutrinoless beta-equilibrium 
(or beta-stable) state of cold matter in the crust below 
the very thin envelope of a neutron star). The typical elec- 
tron fraction of the ejecta is found to be around 0.02-0.04 
(Fig. |lj), much lower than determined by Freiburghaus 
et al. (1999) as favorable for an r-processing that leads to 
a solar-like abundance distribution. We remind the reader 
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Fig. 15. Maximum temperature on the grid as a function 
of time for Models B64 and B64c 



here, however, that the original EoS by Lattimer & Swesty 
(1991) does not provide values below Y^ = 0.03, and we 
use a simple extrapolation of the EoS in this regime. In 
fact, the exact value of Ye in the neutron star crust is not 
well known because of uncertainties in the physics, for ex- 
ample associated with the nucleon symmetry energy or 
with phase transitions due to isospin states in case of very 
low Ye- In addition, the current grid zoning of our three- 
dimensional models does not allow us to resolve the thin 
crust and envelope of the cold initial neutron stars, where 
the state of neutrinoless beta-equilibrium corresponds to 
a positive gradient of Y^ and thus higher values of Yg to- 
wards the surface. 

The question of dynamical mass loss is not unequiv- 
ocally answered on grounds of our simulations or those 
of Rosswog et al. (1999, 2000). A source of uncertainty 
is the nuclear EoS which determines the structure and 
properties of the merging neutron stars and the dynam- 
ics of the merging process. Rosswog et al. (2000) found 
that the mass loss is sensitive to the stiffness of the EoS. 
In particular the conditions in the supranuclear regime 
are highly uncertain. EoSs which have been developed for 
supernova simulations and are particularly suitable for hy- 
drodynamic modeling as described here (because they pro- 
vide all the required information in essentially all interest- 
ing regions of the parameter space), in particular the EoS 
of Lattimer & Swesty (1991) and Shen et al. (1998), do 
actually not include a detailed microphysical description 
of the regime beyond about twice the density of the nu- 
clear phase transition. New hadronic degrees of freedom 
(e.g., a phase with pions, kaons or hyperons) or quark 
matter could be present there and would soften the EoS. 
While probably not crucial for stellar core-collapse and 
supernova simulations, this density regime determines the 
cooling evolution of new-born neutron stars and should 
also affect the merging of neutron stars. For example, if 
the supranuclear EoS is sufficiently soft, the maximum 
mass of stable (nonrotating) neutron stars can be as low 
as 1.5-1.6 Mq. This possibility cannot be excluded, neither 
on grounds of theoretical calculations of the state of mat- 



22 



M. Ruffert and H.-Th. Janka; Coalescing neutron stars - a step towards physical models 



ter at supranuclear densities (e.g., Weber 1999, Heiselberg 
& Pandharipande 2000), nor on grounds of observed neu- 
tron star masses. Although rapid and differential rotation 
can have a significant stabilizing influence (Baumgarte et 
al. 2000), such a soft EoS would probably not allow the 
merger remnant to escape the collapse to a black hole on 
a dynamical timescale (Shibata & Uryii 2001). Therefore 
the question will have to be investigated in more detail by 
future general relativistic simulations, whether in this case 
mass ejection, which occurs with some time delay after the 
final plunge, can still take place. 

4.3. Thermal evolution and neutrino emission 

Computing the temperature evolution is crucial for calcu- 
lating the neutrino emission of the merging neutron stars. 
Unfortunately, there are a number of numerical problems 
which affect the accuracy of the temperature determina- 
tion. For our code and simulations, these problems have 
already been addressed in Sects. || and ^. 

When doing simulations with a finite difference 
scheme, the need to assume a medium with a finite density 
on the computational grid constitutes a problem. In a non- 
rotating frame of reference (or a frame of reference which 
does not move with the same angular velocity as the stars), 
the neutron stars move through this dilute medium with 
a high relative velocity. This leads to shock heating at the 
forward surfaces of the neutron stars. We reset the tem- 
perature to lower values after each time step for densities 
below a certain threshold value so that only dilute mate- 
rial is affected by this procedure. We also detect artificially 
heated, isolated grid zones and reduce the local tempera- 
ture by averaging over neighbouring, well behaved zones. 
Since only relatively few cells with a rather low density are 
involved, this procedure does not introduce a noticable vi- 
olation of the conservation of the total energy during the 
simulations (Fig. Moreover, the affected volume and 
mass decrease when the resolution is improved. 

Another, more severe problem is connected with the 
fact that cool neutron star matter is very degenerate, i.e., 
the thermal energy contributes only a minor fraction to 
the internal energy. Since our hydrodynamics code com- 
putes the specific internal energy, from which the tempera- 
ture is determined, as the difference between internal plus 
kinetic energy of the fluid and its kinetic energy, any small 
error introduced there leads to sizable perturbations in the 
temperature. The thus calculated temperature cannot be 
very accurate. For this reason, we decided to repeat some 
of our simulations by using an additional entropy equa- 
tion, which allows us to follow the thermal history without 
the described sources of noise (see Sect. || for more details) . 
Note that the entropy equation serves only for computing 
the temperature, but does not replace the energy equa- 
tion as one of the conservation laws which are solved for 
describing the hydrodynamical flow. The neutrino source 
terms in the latter equation, however, are computed with 
the temperatures as obtained from the entropy equation. 



Models B32 and B32w in Tables and | are two 
cases which were computed with the same grid resolution, 
but in the second model the additional entropy equation 
was used. Both simulations were started with nearly the 
same temperatures in the neutron stars. The additional 
Model B32c represents a case where the initial central tem- 
perature was chosen to be as low as 0.05 MeV instead of 
6.76 MeV in Model B32 or 5.44 MeV in Model B32w. The 
values for the different entries in the tables show that the 
dynamic quantities and the gravitational-wave emission 
are nearly the same. The quantity that is most sensitive 
to the different treatments is the mass which can be- 
come unbound from the system. Despite of the different 
thermal evolutions, even the overall properties of the neu- 
trino emission at the end of the computations are rather 
similar. 

For discussing more details, let us now consider two 
better resolved calculations. Models B64 and B64c. These 
two exemplary models differ by the use of the entropy 
equation in combination with a low initial central tem- 
perature in Model B64c, namely 0.05 MeV instead of 
6.25 MeV in the neutron stars of Model B64 (Table p. 
The lower initial temperature is more realistic than our 
usually "warm" initial conditions, because viscous heating 
prior to the merging is unlikely to achieve temperatures in 
excess of 10^ K (Bildsten & Cutler 1992, Kochanek 1992, 
Lai 1994). Starting with a cold initial state requires the 
use of the entropy equation (see above). 

Figure ^ shows the maximum temperatures on the 
grid as functions of time for Models B64 and B64c. Indeed, 
the maximum temperatures evolve significantly differ- 
ently. This is confirmed by the upper panels of Fig. ^ 
which show the temperature in the equatorial plane of 
Model B64c for two times, which can directly be com- 
pared with snapshots given for Model B64 in Fig. |. Note 
that in Model B64 the maximum temperature during the 
early phase is limited by the spike-correction procedure, 
whereas such a temperature reduction is not applied to 
the "entropy-temperature" ' in Model B64c. 

The different temperatures in both simulations affect 
the neutrino source terms for lepton number and energy, 
which are very sensitive to the gas temperature (see the 
Appendices in Ruffert et al. 1996). Since the neutrino 
source terms enter the hydrodynamics equations and the 
continuity equation for the electron lepton number, they 
could in principle have consequences for the dynamical 
evolution of the models. In reality, however, the neutrino 
emission is irrelevant for the dynamics, because the asso- 
ciated energy is too small on the short timescale of the 
calculation. 

Although Model B64c has higher peak temperatures 
during the early phase of the evolution, the average 
temperatures are considerably lower than in Model B64 
(Figs. H and [l6|) . These lower temperatures reduce the 
(anyway small) changes of the lepton number due to the 
neutrino emission in the expanding tidal tails. The pan- 
els in the second row of Fig. |l^ confirm that the lepton 
number remains essentially unchanged from its value of 
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Fig. 16. Temperature and electron fraction in the orbital plane of Model B64c for two different times. In this model 
the entropy equation was used to follow the temperature evolution. The temperature contours are spaced linearly with 
levels at 1, 2, 4, 6, 8 and 10 MeV, the Yg contours are given in steps of 0.02, starting with a value of 0.02. The plots 
should be compared with the corresponding panels in Figs. ^ and ^ 



about 0.02 near the surfaces of the original neutron stars 
(compare also with Fig. 

The temperature distributions in the orbital plane of 
Models B64 and B64c at the end of the simulations (see 
Figs. Hand 16) reveal that only moderate differences occur 



in the cloud of low-density material [p ^ 10^^-10"'^^ gcm~^ 
at r ^ 50 km) which surrounds the much denser and very 
compact core of the merger remnant. Shock heating has in- 
creased the temperatures in this cloud to values between 
1 MeV and 5-6 MeV in both models. Much more pro- 
nounced differences between the models can be found in 
the interior of the compact core, where Model B64 is sig- 
nificantly hotter [T > 10-15 MeV), mainly due to the dis- 
sipation of kinetic energy by numerical viscosity. The cor- 



responding heating depends on the shear motions, which 
are particularly strong during the phase when the neutron 
stars merge, and later during the numerous revolutions of 
the rapidly spinning central core within the surrounding 
layer of gas. 

Most of the neutrino emission comes from the outer 
regions of dilute gas, because the dense core of the rem- 
nant is much less transparent to neutrinos. For this reason, 
the luminosities and mean energies of the emitted neutri- 
nos and antineutrinos of all flavors are rather similar for 
Models B64 and B64c (Fig. |l^ , in particular towards the 
end of the simulations, when the extended gas cloud has 
reached a quasi-steady state. Before the tidal arms ex- 
pand and shock heating and viscous shear has raised the 
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Fig. 17. Mean energies (left) and luminosities (right) of Vf, (label "e"), Ve (label "a") and heavy-lcpton neutrinos and 
antineutrinos, Vx (label "x"), as functions of time for Models B64 (top) and B64c (bottom). The luminosity of Vx 
includes contributions from muon and tau neutrinos and antineutrinos 



temperatures in the outer parts of the merged stars, the 
neutrino emission of the initially cold Model B64c in fact 
stays on a relatively low level (~ 10^^ ergs~^). During this 
phase the mean energies of the emitted neutrinos show 
large fluctuations. 

The similarity of the neutrino production in the late 
phase of Models B64 and B64c is also visible in Fig. |l8| , 
where the total loss rates of neutrino energy per unit vol- 
ume are displayed for both models. One can see that the 
dominant contribution to the emission of neutrinos of all 
flavors stems from the outer parts of the merger remnant 
and not from the very dense core. In both models, the 
location of the neutrinospheres, defined where the opti- 
cal depth perpendicular to the orbital plane (i.e., effec- 
tively the minimum in all three coordinate directions of 
the cartesian grid) drops below unity (Eq. 7 in Ruffert & 
Janka 1999), is also very similar. In the orbital plane the 
neutrinosphere radius is around 60-70 km. Perpendicular 
to the equator the density gradient is very steep and the 
neutrinospheres reach up to a vertical distance of about 



20 km. The neutrinospheres of electron neutrinos {ve), 
electron antineutrinos (pe) and of the heavy-lepton neu- 
trinos (/i and T neutrinos and antineutrinos, Vxi are pro- 
duced with the same rates and "see" essentially the same 
opacity) are very close to each other. The neutrinosphere 
of Vx is slightly smaller because there is no contribution 
of charged-current neutrino-nucleon interactions to the Vx 
opacity. 

The typical total neutrino luminosities at the end of 
the computed evolution are of the order of a few times 
lO^^ergs"^ (Fig. |l^). Within about 10ms the models ra- 
diate an energy equivalent of roughly 10~^ Mq in neutri- 
nos, which is more than an order of magnitude less than 
in gravitational waves (compare Figs. |l9| and |l^). These 
neutrino luminosities and energies are considerably larger 
than those given by Ruffert et al. (1997a) for compara- 
ble models. The reason for this result is the use of the 
much larger grid in the current simulations. The previ- 
ous models suffered from the problem that a significant 
amount of matter left the computational volume when the 
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Fig. 18. Total neutrino energy loss rates (in ergcm^'^s^^) from the matter in the orbital plane and perpendicular 
to the orbital plane for Model B64 (top) and Model B64c (bottom) at the end of the computed evolution. The plots 
show the central region of the computational grid with the neutrinospheres of v^i i^e, and i>e marked by dashed lines 
(in this order from inside outward). The contours are spaced logarithmically in steps of 0.5 dex 



tidal tails expanded. In the present calculations this gas is 
wrapped up to form the shock-heated envelope around the 
dense core. Note that the position of the neutrinospheres 
is therefore outside of the computational volume that was 
employed by Ruffert et al. (1997a). 

The largest contribution to the neutrino emission 
comes from electron antineutrinos, because positron cap- 
tures on free neutrons dominate electron captures on pro- 
tons. In the decompressed and hot neutron star matter, 
which forms the envelope around the merger core, neu- 
trons are very abundant and the electron degeneracy is 
moderate, for which reason electron-positron pairs are 
present in large numbers. Because /Lt and r neutrinos and 
antineutrinos are produced only by pair reactions (mainly 



+e+ — !■ i^ + f') — but not by electron/position captures 
on protons/neutrons, which are the dominant processes 
for generating v,, and Pe, respectively — their emission is 
extremely sensitive to the temperature (the energy pro- 
duction rate scales with T^). Only in very hot models do 
they contribute to the neutrino emission on a significant 
level. For example, the difference of the total neutrino lu- 
minosities of Models B64 and B64c is almost entirely due 
to the emission of heavy- lepton neutrinos; the final lu- 
minosities of electron neutrinos and antineutrinos are es- 



sentially the same in both models (Fig. 19). The typical 
mean energies are around 11-13 MeV for the radiated Ve, 
15-19 MeV for 9^, and 22-27 MeV for (Table |). 
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Fig. 19. Total neutrino luminosity (summed up for all flavors) and cumulative energy emitted in neutrinos as functions 
of time for different models 



Figure ^ gives the energy loss rates per unit volume in 
the orbital plane of Model B64 at the end of the simulated 
evolution for all neutrino types (vg, f'e, and Vx) individu- 
ally and for the sum of all contributions. The orbital plane 
is shown out to the boundaries of the computational grid. 
Although the three neutrinospheres nearly coincide, one 
notices clear differences of the spatial distribution of re- 
gions where the different types of neutrinos are produced. 
While Ue are rather uniformly emitted from most of the 
matter, there are definite hot spots which radiate large 
amounts of Ve and Vx- This emphasizes the temperature 
sensitivity of the corresponding neutrino production pro- 
cesses, which require the presence of positrons in the stel- 
lar medium. 

We evaluated two of our models, A64 and B64, for the 
energy deposition by the annihilation of neutrinos and an- 
tineutrinos to electron-positron pairs, v + D ^ +e'^ , as 
described by Ruffert et al. (f997a) and Ruffert & Janka 
(1999). The results at the end of the computed evolution 
(where the neutrino luminosities are highest) are plotted 
in Fig. |2^. The largest energy deposition rates per unit 
volume are found above the poles of the merger remnant 
with peak values in excess of 10^^ ergcm^^s^^. The num- 
bers for the total rate of energy deposition in the vol- 
ume where the gas density is below 10^^ gcm^^ are about 
9 X 10^^ ergs"^ in case of Model A64 and 7 x lO'^^ ergs"^ 
for Model B64 (Table |). These rates are more than 20- 
30 times bigger than calculated by Ruffert et al. (1997a). 
This is mainly explained by the much higher neutrino lu- 
minosities in our current simulations, which enter the an- 
nihilation rate quadratically. (To some degree the effect 
is also caused by the different geometry, compare Fig. ^ 
with Fig. 16 in Ruffert et al. 1997a.) 

The largest and by far dominant part of the energy 
deposition occurs in the polar regions. Here, however, also 
the energy loss rates by neutrino emission reach max- 
ima (see the panels in the right column of Fig. |l^), be- 



cause high temperatures are present in a region where the 
gas is still dense, but the density gradient is very steep 
(Fig. 1^). Therefore neutrinos are generated in large num- 
bers and can easily escape to the transparent regime, as 
suggested by the fact that the neutrinospheres cross the re- 
gions of peak emission above the poles of the compact core 
(Fig. |l^). With values up to more than 10'^'^ ergcm~'^s~^, 
the rate of energy loss exceeds the energy input rate by 
annihilation by more than a factor of about 10. Therefore 
the energy which is transferred to the stellar plasma will 
efficiently and immediately be reradiated by neutrino pro- 
duction processes, and the numbers for the total energy 
deposition rate in Table ^ overestimate the net heating ef- 
fect by a large factor. The conclusion drawn by Ruffert et 
al. (1997a) and Janka & Ruffert (1996) remains valid: The 
energy deposited by neutrino-antineutrino annihilation in 
the neutrino-transparent, low-density plasma before, dur- 
ing, and immediately after the merging of two neutron 
stars is not sufficient to explain the energetics of typical 
gamma-ray bursts. 

4.4. Black hole formation and accretion 

While the emission of gravitational waves peaks right at 
the moment when the two neutron stars fall into each 
other, the neutrino luminosity does not reach a very high 
level before the tidal tails have been inflated and wrapped 
up to the hot gas cloud that finally surrounds the dense 
core. Even then, neutrino-antineutrino annihilation is not 
an efficient mechanism to provide the energy for a gamma- 
ray burst, because electron and positron captures on free 
nucleons extract the deposited energy extremely rapidly 
from the dense gas that is present also above the poles of 
the merger remnant. 

The situation changes, when the core of the remnant 
collapses to a black hole. Our basically Newtonian sim- 
ulations, however, do not yield evidence whether and if 
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Fig. 20. Neutrino energy loss rates (in ergcm~'^s~^) from the matter in the orbital plane of Model B64 at the end of 
the computed evolution. The rates for v,, (top left), D,, (top right), the sum of muon and tau neutrinos and antineutrinos 
(bottom left) and for neutrinos and antineutrinos of all flavors are given logarithmically, with steps of 0.5 dex. The 
dashed lines in the last plot indicate the positions of the neutrinospheres in the orbital plane 



so, when such a collapse occurs. Assuming that it hap- 
pens, Ruffert & Janka (1999) continued the simulation of 
Model B64 for several ms to investigate the effects on the 
surrounding gas cloud. The black hole was represented by 
a gravitating "vacuum sphere" , and the time of the grav- 
itational instability was treated as a free parameter. The 
results for the dynamical evolution and the neutrino emis- 
sion, however, turned out to be rather similar, indepen- 
dent of whether the black hole was assumed to form only 
~ 2 ms or as late as 10 ms after the start of Model B64. 

A comparison of Fig. 15 in Ruffert & Janka (1999) 
with Fig. in the present paper reveals that the baryonic 
matter above the poles of the compact core falls into the 
newly formed black hole very quickly. Within milliseconds. 



a funnel along the rotation axis is cleaned from baryons, 
and the density decreases to a value near our numerical 
lower limit of the density. With most of the energy depo- 
sition by ncutrino-antineutrino annihilation taking place 
in that region, this appears to be a favorable situation for 
a baryon-poor, potentially relativistic outflow to develop, 
which might later on produce a gamma-ray burst through 
dissipation of the mechanical kinetic energy into radiation. 

Besides the formation of the black hole, such a scenario 
requires that a significant amount of matter remains in a 
disk around the black hole. The accretion of this matter 
on a timescale much longer than the dynamical timescale 
allows for a high efficiency of the conversion of rest-mass 
energy of the accreted gas to neutrinos. Typically several 
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Fig. 21. Energy deposition rate (in crgs^^cm^'') by the annihilation of neutrinos and antineutrinos to electron- 
positron pairs in the surroundings of the merger remnant for Models A64 (left) and B64 at the end of the computed 
evolution. Azimuthally averaged rates are shown in a plane perpendicular to the orbital plane. The corresponding solid 
contours are spaced logarithmically in steps of 0.5 dex. The dotted lines indicate the (azimuthally averaged) isodensity 
levels, also spaced logarithmically with steps of 0.5 dex. The rates are evaluated only in regions where the density is 
less than 10^^ gcm^"^ 



per cent efficiency in case of a Schwarzschild black hole 
and up to several ten per cent are possible for a Kerr 
black hole which accretes matter from a corotating (thin) 
disk (Thorne 1974, Shapiro & Teukolsky 1983, Popham et 
al. 1999, Li & Paczynski 2000). Even more energy is avail- 
able when the rotational energy of the black hole is tapped 
by means of magnetic fields (Blandford & Znajek 1977, 
Li 2000, Lee et al. 2000). 

If a black hole forms from most of the mass of the 
merger remnant, our hydrodynamical models allow us to 
estimate the mass which ends up in an accretion disk 
around this black hole. Assuming that the disk is sup- 
ported mainly by centrifugal forces, we use the criterion 
that the specific angular momentum of the gas should be 
larger than the Keplerian angular momentum at three 
Schwarzschild radii, i.e., at the location of the inner- 
most stable circular orbit for a nonrotating black hole: 
j > VGGAf/c, where M is taken to be the total mass on 
the grid. The corresponding gas masses, Md, at the end of 
our simulations are listed in Table |^, and are indicated for 
Models B64 and D64 in Fig. |22[ Typically, several hun- 
dredths of a solar mass up to a few tenths of a solar mass 
fulfill this condition. The largest numbers are obtained for 
corotating models (B, S, D), where the gas has the high- 
est specific angular momentum. The continuation of the 
simulation of Model B64 by Ruffert and Janka (1999) con- 
firmed that these estimates are in reasonably good agree- 
ment with the gas mass which finally orbits around the 
black hole when a quasi-stationary state is reached. 

For the irrotational A-case (and in the cases where the 
neutron stars started out with opposite spin directions or 
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Fig. 22. Masses on the grid with a density below 
101° gcm-3 ("Mio"), below 10" gcm-^ ("Mn"), and 
below 10^^ gcm-3 ("M12"), respectively, for Models B64 
(solid lines) and D64 (dashed lines) as functions of time. 
The thin curves mark the estimated disk masses, Md, 
which encompass all gas with a specific angular momen- 
tum larger than the Keplerian angular momentum at three 
Schwarzschild radii of the merger remnant (i.e., of the to- 
tal mass on the grid) 

spins in counterrotation to the orbit), only a few per cent 
of the total rest mass of the binary system can gather 
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enough angular momentum by hydrodynamical processes 
to be able to stay in a disk around a black hole at the 
end of our ~ 10 ms of computed evolution. These esti- 
mates are confirmed by three-dimensional simulations of 
neutron star mergers in full general relativity (Shibata & 
Uryu 2000, 2001). Of course, if the collapse to a black hole 
occurs immediately after the merging of the two neutron 
stars, before the tidal tails and the extended envelope of 
the merger remnant have a chance to form, there is no 
time for any transport of angular momentum by hydrody- 
namic interaction, and very little gas, if any at all, will be 
able to remain in a disk. 

On grounds of our current hydrodynamical simula- 
tions, solving the time-dependent Euler equations, we can- 
not draw conclusions on the further development of the 
accretion "disk" or, better, of the extended torus. Once 
the gas has settled into a quasi-steady state around the 
newly-formed black hole, its later destiny will be driven 
by the radial transport of angular momentum, which is 
mediated by viscous forces, and the torus structure and 
internal conditions will be determined by the energy (and 
lepton) number loss through neutrino emission. Magnetic 
fields can play an important role, too. In our simulations, 
without taking into account the effects of the physical disk 
viscosity (the numerical value of which cannot be deter- 
mined from first principles and thus would have to be con- 
sidered as a free parameter within the Navier-Stokes equa- 
tions), the subsequent evolution is governed by the action 
of the numerical viscosity, which is not under our direct 
control, but can be varied only indirectly by changing the 
grid resolution. In addition, numerical viscosity does not 
have the same properties as the physical viscosity, e.g., it 
does not necessarily conserve the angular momentum. 

Considering the situation at the end of our models we 
therefore could only speculate about the long-time evolu- 
tion of the accretion torus (Ruffert & Janka 1999, Janka et 
al. 1999). With a given value for the mass accretion rate by 
the black hole, for example, we obtained an estimate of the 
torus lifetime. Moreover, the mass, density, and tempera- 
ture of the torus define its neutrino emission properties, 
and the calculated neutrino luminosity (extrapolated for 
the estimated accretion time) allows one to come up with 
numbers for the efficiencies of conversion of gravitational 
potential energy to neutrinos and to the electron-positron 
pair plasma by neutrino-antineutrino annihilation. 

Typical temperatures within the torus are a few MeV, 
the maximum temperatures between about 5 MeV and 
roughly 10 MeV. A significant amount, if not most of the 
torus mass has a density above 10^^ gcm~^ (Fig- ^ and 
Table "Massive" tori, i.e., those with a mass of more 
than 0.1 Mq, are optically thick to neutrinos, whereas 
tori with masses of only a few 0.01 M© are close to neu- 
trino transparency (see Fig. 30 in Ruffert & Janka 1999). 

Whether the dense core of the merger remnant col- 
lapses to a black hole, and if so, on what timescale this 
happens, is a complex question, which requires not only a 
general relativistic treatment, but depends on a number of 
additional aspects, e.g., the mass and compactness of the 
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Fig. 23. Azimuthally averaged angular velocity in the or- 
bital plane as a function of the distance from the system 
axis. The plot shows the inner regions of the merger rem- 
nants of Models A64, B64, 064, and D64 at the end of the 
computed evolution 

initial neutron stars, the properties of the nuclear EoS, 
and the rotation of the post-merging object. Shibata & 
Uryu (2000, 2001), by performing three-dimensional dy- 
namical simulations in full general relativity, find that the 
product of a neutron star merger is sensitive to the ini- 
tial compactness of the neutron stars, which is defined 
as the ratio of the Schwarzschild radius of the neutron 
star to its actual radius. This quantity increases when the 
mass of the neutron star approaches the limiting mass of 
a spherical star in isolation. For sufhciently compact stars 
a black hole is formed on a dynamical timescale, as the 
compactness descreases, the formation timescale becomes 
longer and longer. (The corresponding critical mass of the 
binary relative to the mass limit of a nonrotating, single 
neutron star varies with the compressibility of the EoS.) 
For less compact differentially rotating "supra- 

massive" neutron star (Cook et al. 1992, 1994a,b) forms, 
which first has to become a rigidly rotating body (on a sec- 
ular timescale by viscous dissipation) or has to lose some of 
its angular momentum by dynamical processes (e.g., mass 
stripping, bar-mode instability) or secular processes (e.g., 
gravitational waves, mass loss by winds, magnetic fields), 
before the gravitational instability can set in (Shibata et 
al. 2000). 

Thus the evolution leads to a black hole very quickly 
only if the total rest mass of the binary system is suffi- 
ciently larger than the maximum rest mass of a spherical 
star in isolation. The exact factor depends on the stiffness 
of the EoS as well as on the angular momentum of the 
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binary system. The latter determines the rotation of the 
core of the merger remnant, which can remain stable al- 
though it has a mass that is significantly larger than the 
mass limit of spherical or rigidly rotating neutron stars 
because of the supporting effect of rapid, differential rota- 
tion (Baumgarte et al. 2000). 

Our models have a dimensionless rotation parameter 
a = Jc/{GAi'^) (J is the total angular momentum, M 
the total rest mass of the system), which initially is be- 
tween 0.64 for Model C64 and 0.98 for Model S64. During 
the evolution, the mass changes only slightly, because lit- 
tle gas escapes from the system. In contrast, the angular 
momentum decreases significantly due to the emission of 
gravitational waves (and partly due to angular momentum 
which is carried away by the ejected gas). The final values 
of the parameter a are therefore lower, between about 0.5 
for Model C64 and 0.75 for Model S64 (see the columns oi 
and Of in Table 2 of Janka et al. 1999; the actual numbers 
might be up to ~ 10% larger, because one should apply 
corrections for the numerical loss of angular momentum) . 
This suggests that the merger remnant has an angular 
momentum below the critical limit that is possible for a 
Kerr black hole {a — 1). Of course, a definite statement of 
this kind is not possible without a relativistic treatment, 
which considers M to be the gravitational mass instead 
of the rest mass, and includes the mass reducing effect of 
the energy loss by gravitational waves. The latter, how- 
ever, is small. According to the quadrupole formula less 
than ~ 1% of the mass of the system is emitted in gravi- 
tational waves (Fig. |ll|). Therefore the trend seen in our 
Newtonian calculations should be right, and the subcriti- 
cal rotation of the merger remnant should even be quanti- 
tatively correct if one does not start with an initial value 
of a that is much larger than unity. 

It is worth a final remark that the compact core of 
the merger remnant at its "surface" near 25 km rotates 
with an angular velocity that is only slightly lower than 
the Keplerian frequency, J^k = a/ GM/r^ ~ (4. ..5) x 
10^ s~^. The degree of differential rotation, however, varies 
strongly between the different models and depends on the 
initial neutron star spins and the mass ratio of the neu- 
tron stars. Figure |2^ shows the azimuthally averaged an- 
gular velocities in the orbital plane for the core of the 
remnant in different models at the end of the computed 
evolution. Whereas in the irrotational, symmetric case. 
Model A64, the core rotates moderately differentially with 
radially decreasing mean angular velocity, the asymmetric 
Model D64 reveals the opposite trend, and the corotating 
case. Model B64, is near rigid rotation. The most dra- 
matic differential rotation is found in the inner region of 
Model 064, where the mean angular velocity drops from 
roughly 14 x 10^ s^^ near the center to a value around 
3.5 X 10'^ s~^ between 15 km and 30 km distance from 
the system axis. We would like to emphasize, however, 
that except for Model 064, the different grid zones at a 
fixed equatorial radius show large fluctuations of the corre- 
sponding value of f2, because the merger remnants are still 
significantly perturbed and perform violent oscillations. 



5. Summary, discussion, and conclusions 

We have presented results of neutron star merger simula- 
tions that were performed with a new version of our nu- 
merical code, which was significantly improved compared 
to the original one (Ruffert et al. 1996, 1997a), mainly by 
the introduction of several levels of nested grids. These 
allow for a better resolution of the stars near the grid cen- 
ter on the one hand, and a larger computational volume 
on the other. Besides recomputing our previously pub- 
lished models with the new code, we varied the neutron 
star parameters (masses, mass ratios, spins) and computed 
models with different maximum resolution. The lower res- 
olution on the finest grid was similar to the grid in our 
older calculations, whereas the current "standard mod- 
els" have a factor of two higher resolution in all three 
cartesian directions. Moreover, we tested the accuracy (or 
uncertainty) of the temperature calculation by replacing 
the energy (i.e., the sum of kinetic and internal energies 
in our code) by the entropy as the basic variable to follow 
the temperature evolution in time. 

Using the energy for deriving temperatures has the 
disadvantage that the thermal energy is only a small con- 
tribution to the internal energy. The latter is dominated 
by the degeneracy energy of the fermions. This leads to er- 
rors in the temperature determination when the internal 
energy is not very accurate. With the entropy equation 
one also reduces, although cannot eliminate, the problem 
that unphysical shocks at the interface between the neu- 
tron stars and the surrounding medium occur and cause 
an overestimation of the temperature, in particular be- 
fore the merging. Even more, the entropy equation gives 
one control over the effects of shear and bulk viscosity 
on the thermal evolution, whereas the temperature calcu- 
lated from the energy equation does not allow this direct 
control, because the dissipative effects of numerical vis- 
cosity in the code can be influenced only by changing the 
grid resolution. Of course, solving the entropy equation 
as a supplementary equation in addition to the conserva- 
tion laws of mass, momentum, energy, and lepton num- 
ber, which still describe the evolution of the stellar fluid 
(and where the neutrino source terms were evaluated with 
the temperature as obtained from the entropy equation), 
is not a hydrodynamically fully consistent approach. We 
therefore do not consider this procedure as the necessar- 
ily more accurate calculation, but as an attempt to test 
the sensitivity of our results to effects that are associated 
with the uncertainties of the temperature determination 
discussed above. It is reassuring that despite of significant 
differences of the thermal evolution our main results show 
rather little variation. 

In fact, this work was also strongly motivated by 
the wish to investigate and outline (at least some) ma- 
jor uncertainties of (not only our) current neutron star 
merger simulations. Such uncertainties have a bearing on 
the possibility to draw model-based conclusions on the 
gravitational-wave emission, a potential connection with 
gamma-ray bursts, and the implications for the production 
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of heavy elements in our Galaxy. The models presented in 
this paper improve our previous calculations with respect 
to numerical resolution and reduced mass loss through 
the outer boundaries of the significantly enlarged compu- 
tational grid. They are intended to serve for comparisons 
with future general relativistic simulations. 

5.1. Main results 

In summary, our main results and their implications are 
the following. The details of the gravitational-wave sig- 
nal, i.e., the primary and subsequent maxima of the lu- 
minosity, the total radiated energy, and the structure of 
the wave amplitude in particular during the post-merging 
phase, when the core of the merger remnant is in a highly 
perturbed state and performs violent oscillations, exhibit 
some change with the resolution on the finest grid. A cell 
size of 0.64 km or smaller, corresponding to at least 50 
zones on the diameter of the initial neutron star, seems 
desirable. With less resolution, the effects of numerical vis- 
cosity and the associated loss of angular momentum grow 
to an unacceptable level and characteristic features of 
the gravitational-wave emission are noticably influenced. 
This limits the possibility to deduce important informa- 
tion from a possible future wave measurement. The grav- 
itational waves which accompany the final plunge of the 
neutron stars, carry information about the masses, com- 
pactness and the spins of the neutron stars and thus allow 
for conclusions on the nuclear EoS. The post-merging sig- 
nal yields evidence about the destiny of the merger rem- 
nant and can also be used to extract information about 
the internal state of neutron stars. 

While the gravitational-wave emission is strongest 

when the two neutron stars plunge into each other, the 
neutrino emission rises only gradually afterwards. It ap- 
proaches a saturation level towards the end of our simu- 
lations, when the tidal tails have been wrapped up to a 
cloud of shock- (and shear-)heated gas that surrounds the 
compact and much denser core of the merger remnant. 
Because of the more extended computational grid, which 
is necessary to follow the development of this gas cloud, 
the neutrino luminosities are found to be a factor of 2-4 
higher than in our previously published models (Ruffert 
et al. 1997a). Correspondingly, the energy deposition rates 
by neutrino- antineutrino annihilation in the dilute outer 
layers {p < 10^"'^ gcm~^) of the post- merging object arc up 
to a factor of 20-30 larger. By far the dominant part of this 
energy, however, is deposited in the polar regions, where 
the temperature is high and the scattering depth to neu- 
trinos is rather low. Therefore this energy is immediately 
r eradiated through neutrinos produced by electron and 
positron captures on nucleons. Such an energy transfer to 
a region with large baryon density is therefore inefficient in 
powering a gamma-ray burst, and should drive a baryonic 
wind rather than a relativistic outflow of a baryon-poor 
pair-photon plasma. This presumably high-entropy wind 
(i.e., the medium has low density but comparatively high 



temperature), which expands into vacuum, may have very 
interesting, so far not investigated, implications for ob- 
servable radiation from neutron star merging events, and 
for the enrichment of our Galaxy with heavy elements 
(a more detailed discussion can be found in Ruffert & 
Janka 1999 and Janka & Ruffert 2001). 

The dynamical mass ejection from the merging binary 
varies with the initial spins of the neutron stars. It is 
largest (~ (2. ..4) x 10~^ Mq or roughly 1% of the system 
mass) for corotating systems (a case which is not likely 
to be realized because of the small viscosity of neutron 
star matter, which cannot bring the system into a tidally 
locked state prior to merging; Kochanek 1992, Bildsten & 
Cutler 1992) and smallest (;$ IO^'^Mq) when the stars 
initially counterrotate with the orbit. Only in the for- 
mer case very prominent tidal arms develop through the 
outer Lagrange points and expand away from the center 
of the merger. The use of the larger computational vol- 
ume helped to signiflcantly improve our estimates for the 
amount of mass which can potentially become unbound. 

It is an important aim of numerical models to de- 
termine the thermodynamical conditions and the nuclear 
composition of these ejecta, and their subsequent evolu- 
tion. This will help answer the question whether and how 
r-processing can take place in this matter (Lattimer & 
Schramm 1974, 1976; Symbalisty & Schramm 1982; Meyer 
1989; Davies et al. 1994; Freiburghaus et al. 1999) and 
whether it has contributed to the heavy-element content 
of our Galaxy at a significant level. 

We find that the matter, which is ejected from the tips 
of the expanding tidal tails, stays cool, because it is nei- 
ther heated by shocks nor by viscous friction. In fact, it 
is a problem in hydrodynamical simulations to accurately 
trace the thermal history of the initially cold neutron mat- 
ter (the viscosity is not only too small to achieve tidal 
locking, it is also too small to heat the neutron stars be- 
yond ~ 10® K before merging; Lai 1994). Besides the high 
degeneracy of the medium, numerical viscosity, which is 
present to some (but actually different) degree in all nu- 
merical codes and depends on the resolution, causes prob- 
lems for an accurate calculation of the temperature. In ad- 
dition, the limited numerical resolution does not allow one 
to properly represent the extremely steep density gradi- 
ent near the neutron star surface below a density of about 
lO^^gcm"^. This leads to the necessity of softening the 
density gradient to produce (nearly) hydrostatic condi- 
tions. In our Eulerian, grid-based simulations the neutron 
stars also have to be embedded by a low-density medium. 
When the neutron stars move through this surrounding 
medium, shocks occur at the stellar surfaces in upstream 
direction, which produces artificial heating. Thus locally, 
the temperatures can be overestimated by a large factor, 
although this numerical heating is small compared to the 
maximum temperatures which are reached when the two 
neutron stars plunge into each other. Because of all these 
aspects, some of which arc not specific to a particular code 
but are generic to the physical problem or to hydrody- 
namical simulations with limited resolution, the calculated 
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temperatures are likely to be overestimated and have to 
be interpreted with special caution. Significant progress 
requires a much enlarged numerical resolution. For grid- 
based codes this could be achieved by employing adaptive 
mesh refinement techniques. Another option for improv- 
ing some aspects may be the choice of a rotating instead 
of a fixed frame of reference (New & Tohline 1997, Swesty 
et al. 2000), although the inspiral is so rapid that there 
would quickly be motion of the stars also in a reference 
frame that is initially corotating. 

Our merger simulations followed the variation of Ye 

in response to electron captures by protons and positron 
captures by neutrons, including also the effects of neutrino 
trapping when the stellar medium becomes optically thick 
to neutrinos at sufficiently large densities. Although these 
weak interaction rates are fast at temperatures between 
10^° K and 10^^ K, they are not effective in raising the 
initial electron fraction on the very short timescale of the 
dynamical expansion of the cjccta. Starting with a typical 
equilibrium Ye around 0.02 in the tidal tails, we cannot 
find a growth by more than a factor of 2 3 (to values of at 
most ~ 0.06) before the gas leaves our computational grid. 
At that time the density has decreased to 10^°gcm~"^ or 
less, and the matter has cooled down to (2-3) x 10^° K. 
Since both the temperature and the density drop quickly, 
and more and more nucleons recombine to alpha particles 
and nuclei, we do not expect a significant effect due to 
electron and positron captures during the subsequent ex- 
pansion. This result was obtained although the tempera- 
ture (and therefore the mass fraction of free nucleons) was 
overestimated in our simulations, implying unrealistically 
fast electron and positron captures on the free protons and 
neutrons. In comparative calculations, using the entropy 
equation for following the thermal history of the medium 
(see above) and starting with lower (and thus more real- 
istic) temperatures, we actually cannot detect any change 
from the initial value of Yg- We emphasize that our point 
here is the fact that neutrino processes seem to be unable 
to alter Ye during the rapid decompression of essentially 
cold neutron star matter. The initial Ye of the neutron 
star crust is therefore preserved, although the exact value 
may be EoS dependent and is therefore uncertain. 

5.2. Elements of a possible r-process site 

Cold material in neutron star crusts consists of neutron- 
saturated, very neutron-rich nuclei that are arranged on 
a lattice and immersed in a gas of neutrons and degener- 
ate electrons (e.g., Weber 1999). This region has a very 
low proton-to-nucleon ratio {Ye ~ 0.02; only in a thin skin 
layer of the neutron star, the outer crust, the electron frac- 
tion rises again). As discussed above, it is unlikely that 
the decompressed and expanding matter in the tidal tails 
is heated by dissipative processes (shocks, viscous shear) 
to temperatures where nuclear statistical equilibrium sets 
in (T ^ 5 X 10^ K). Therefore the memory of the ini- 
tial nuclear composition is not erased during this phase 



of the expansion, and the subsequent changes of the nu- 
clear abundances have to be determinc;(l by following the 
beta-decays of the initial ensemble of heavy nuclei (Hilf 
et al. 1974, Lattimer et al. 1977, Meyer 1989). Since heat- 
ing by beta-decays can raise the temperatures to several 
10^ K without, however, necessarily destroying the mem- 
ory of the initial composition (Meyer 1989, Sumiyoshi et 
al. 1998), a self-consistent coupling of the hydrodynam- 
ical evolution with the effects of nuclear decays, includ- 
ing the decay heating and possible (7,p) and (7,n) reac- 
tions, is essential for reliable predictions of the final nucle- 
osynthetic yields. Calculations fulfilling these conditions 
have not been performed so far, and it remains to be seen 
whether this scenario yields a solar system like distribu- 
tion of r-process elements. 

In contrast, Preiburghaus et al. (1999) assumed that 
the ejected gas starts from nuclear statistical equilibrium, 
i.e. with a temperature above ~ 5 x 10^ K, because they 
used the EoS of Lattimer & Swesty (1991), which was 
actually developed for supernova conditions and does not 
contain the physics required to describe cold neutron stars. 
Combining hydrodynamic results of neutron star merger 
simulations with network calculations, they found that for 
proton-to-nucleon ratios around 0.1 rapid neutron cap- 
tures produce an abundance pattern which fits the ob- 
served solar r-process very well for nuclear masses around 
and above the A « 130 peak. However, they considered 
Yg as a free parameter instead of determining it as a re- 
sult of neutrino emission and absorption reactions in the 
hydrodynamical merger model. Moreover, the feedback of 
beta-decay heating on the hydrodynamic evolution of the 
fluid elements was not taken into account. 

Based on our simulations we come to the conclusion 
that the ejected gas stays cool, does not get heated by 
shocks or viscous shear to tcmpcratiires where nuclear 
statistical equilibrium holds, and retains its very low ini- 
tial proton-to-nucleon ratio. Therefore our simulations 
do not yield the conditions which were determined by 
Preiburghaus et al. (1999) as favorable for producing a 
solar-like r-process abundance pattern in the A ^ 130 
mass range. Our models do not provide evidence that their 
values for the initial temperature and composition are re- 
alised in the ejecta from neutron star mergers. Instead, 
our results seem to favor the kind of scenario discussed by 
Lattimer et al. (1977) and Meyer (1989), who considered 
the decompression of initially cold matter from the inner 
crust of neutron stars, where the composition is dominated 
by extremely neutron-rich nuclei (heavy metals) that can 
be arranged on a lattice and are immersed in a gas of 
neutrons and relativistic electrons. However, it is unclear 
whether the decompression from such an initial state leads 
to the abundance distribution of rapid neutron capture el- 
ements as observed in the solar system. 

Besides beta-decays, electron and positron captures 
and (7,p) and (7, n) reactions, a detailed discussion of the 
nucleosynthesis in the dynamically ejected matter might 
also require taking into account the interaction with the 
intense fluxes of energetic neutrinos from the merger rem- 
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nant. Neutrinos absorbed by nucleons and scattered by 
nuclei in the outflowing gas may heat and accelerate the 
gas, may change the proton-to-neutron ratio and may re- 
process the heavy elements by inducing nuclear transmu- 
tations. Although the neutrino emission from the rem- 
nant rises only gradually after the tidal tails have formed, 
and the timescale of the expansion of the gas away from 
the merger site is very short (of the order of the escape 
timescale, which is roughly tcxp ~ r/vcsc ^ 1ms), the 
number of neutrino-nucleon interactions can be estimated 
to be significant. An order of magnitude evaluation shows 
that about 10% of the nucleons might react with neutri- 
nos: 

/ dr ~ 0.1 a (e,)2o (Afg R^jV^^^ , (H) 

where IZ is the reaction rate per neutrino per unit of 
volume, (e^)2o the mean energy of the emitted neutrinos 
measured in 20 MeV, Li/,53 the total neutrino luminosity 
in lO^'^ergs"^, M3 the mass of the remnant normalized 
to 3M0, and i?i,7 the initial radius in units of 10^ cm. 
The radial velocity was assumed to be roughly given by 
v{r) — y^2GM /r, and the neutrino interaction cross sec- 
tion was approximated by the absorption cross section of 
fe and Ue on nucleons, a ^ 10~^'^(ejy/l MeV)^ cm^, be- 
cause electron neutrinos and antineutrinos dominate the 
neutrino emission of the merger remnant. The quantity a 
is a geometrical factor and is less than unity. It accounts 
for the fact that the neutrinos are radiated mainly per- 
pendicular to the orbital plane (because in this direction 
the scattering depth is smaller), whereas the dynamically 
ejected gas moves away from the system in the plane of 
the binary orbit. 

The amount of dynamically ejected material varies 
strongly with the neutron star spins and is largest for 
the improbable case of corotating systems. Moreover, it 
is sensitive to the stiffness of the unknown nuclear EoS 
(Freiburghaus et al. 1999, Rosswog et al. 2000). Because 
momentum has to be transferred by hydrodynamical pro- 
cesses before gas can be expelled, mass ejection may even 
be suppressed by a quick collapse of the merger remnant 
to a black hole, a possibility which depends on the neutron 
star equation of state, the masses of the merging stars, and 
the angular momentum of the binary system (Shibata & 
Uryii 2001). In view of these fundamental uncertainties, 
current models are unable to yield quantitatively meaning- 
ful numbers for the contribution of neutron star mergers 
to the metal enrichment of our Galaxy, even if the theo- 
retical estimates of the merger rates were reliable (which 
in fact is not the case). 

5.3. Concluding remarks 

We have presented results from state-of-the-art simula- 
tions of neutron star mergings with the most advanced 
version of our hydrodynamics code, and have discussed 
these results and their limitations. Of course, our calcula- 
tions are far from being satisfactory. Relativistic simula- 



tions are necessary to address the question whether a black 
hole forms and if so, on what timescale it happens. This 
is important not only for predictions of the gravitational- 
wave signal and the mass which can become unbound dur- 
ing the merging. It is also important for studying the im- 
plications of neutron star mergers for the origin of heavy 
elements and has consequences for potentially observable 
signals connected with the neutrino- and photon-cooling 
phases of the remnant, which is either a rotating (supra- 
massive), hot neutron star or a black hole which accretes 
matter from a surrounding torus at very high rates. 

So far our treatment of neutrino production and emis- 
sion by using a trapping scheme docs not allow us to study 
the effects of neutrino transport and neutrino energy de- 
position in the outer layers of the merger remnant. The 
latter should drive a baryonic outflow from the massive 
neutron star or from the accretion disk around the black 
hole. Including the feedback of neutrino interactions is es- 
sential for modeling this mass loss and for determining the 
conditions in the wind. 

This discussion shows that the modeling of the merg- 
ing of neutron stars and of the accompanying physical pro- 
cesses is still at its beginning. Current simulations do nei- 
ther fully account for the effects of general relativity, nor 
do they include all the physics which is relevant to come 
up with meaningful predictions of the potentially observ- 
able photon and neutrino emission, or to clarify the role of 
neutron star mergers for the production of heavy elements 
in our Galaxy. Therefore conclusions drawn from current 
hydrodynamical models have to be considered with spe- 
cial reservation. We have outlined a number of aspects 
where improvements and progress in future simulations 
are highly desirable. 

Acknowledgements. It is a pleasure for us to thank Wolfgang 
Keil for transforming Lattimer & Swesty's FORTRAN equa- 
tion of state into a usable table. We would also like to 
thank an anonymous referee for his careful reading and use- 
ful comments. MR is grateful for support by a PPARC 
Advanced Fellowship, HTJ acknowledges support by the 
"Sonderforschungsbereich 375 fiir Astro- Teilchenphysik" der 
Deutschen Forschungsgemeinschaft. The calculations were per- 
formed at the Rechenzentrum Garching on a Cray-YMP 4/64 
and a Cray-EL98 4/256. 

References 

Akmal A., Pandharipande V.R., RavenhaU D.G., 1998, Phys. 

Rev. C 58, 1804 
Ayal S., Piran T., Oechslin R., Davies M.B., Rosswog S., 2001, 

ApJ 550, 846 

Baumgarte T.W., 2001, in Astrophysical Sources of 
Gravit ational Radiati on, ed. J.M. Centrella, AIP Press, in 
press ( br-qc/010104^ ) 

Baumgarte T.W., Shapiro S.L., Shibata M., 2000, ApJ 528, 
L29 

Baumgarte T.W., Cook G.B., Scheel M.A., Shapiro S.L., 
Teukolsky S.A., 1998a, Phys. Rev. D 57, 6181 

Baumgarte T.W., Cook G.B., Scheel M.A., Shapiro S.L., 
Teukolsky S.A., 1998b, Phys. Rev. D 57, 7299 



34 



M. Ruffert and H.-Th. Janka; Coalescing neutron stars - a step towards physical models 



Bildsten L., Cutler C, 1992, ApJ 400, 175 
Blanchet L., Damour T., Schafer C, 1990, MNRAS 242, 289 
Blandford R.D., Znajek R.L., 1977, MNRAS 179, 433 
Bonazzola S., Gourgoulhon E., Marck J. -A., 1999, Phys. Rev. 
Lett. 82, 892 

Bulik T., Belczynski K., Zbijewski W., 1999, MNRAS 309, 629 
Centrella J.M., McMillan S.L.W., 1993, ApJ 416, 719 
Colella R, Woodward RR., 1984, J. Comput. Phys. 54, 174 
Cook G., Shapiro S.L., Teukolsky S.A., 1992, ApJ 398, 203 
Cook C, Shapiro S.L., Teukolsky S.A., 1994a, ApJ 422, 227 
Cook C, Shapiro S.L., Teukolsky S.A., 1994b, ApJ 424, 823 
Cooperstein J., 1988, Physics Reports 163, 94 
Cutler C, Flanagan E.E., 1994, Phys. Rev. D 49, 2658 
Davies M.B., Benz W., Piran T., Thielemann F.K., 1994, 
ApJ 431, 742 

Duncan R.C., Shapiro S.L., Wasserman I., 1986, ApJ 309, 141 
Eichler D., Livio M., Piran T., Schramm D.N., 1989, Nature, 
340, 126 

Faber J. A., Rasio F.A., 2000, Phys. Rev. D 62, 064012 
Faber J.A., Rasio F.A., 2000, Manor J.B., Phys. Rev. D 63, 
044012 

Flanagan E.E., 1999, Phys. Rev. Lett. 82, 1354 
Freiburghaus C, Rosswog S., Thielemann F.-K., 1999, 
ApJ 525, L121 

Fryer C.L., Woosley S.E., Hartmann D., 1999, ApJ 526, 152 
Haensel P., Schaeffer R., 1992, Phys. Rev. D 45, 4708 
Heiselberg H., Pandharipande V., 2000, Ann. Rev. Nucl. Part. 
Sci. 50, 481 

Hilf E.R., Hillebrandt W., Takahashi K., El Eid M.F., Kodama 

T., 1974, Physics Scripta lOA, 132 
Hulse R.A., Taylor J.H., 1975, ApJ 195, L51 
Janka H.-T., Ruffert M., 1996, A&A 307, L33 
Janka H.-T., Ruffert M., 2001, in Proc. Conf. on Stellar 

Collision s and Mergers, ed . M. Shara, ASP Conf. Series, 

in press ( j,stro-ph/0101357 ) 
Janka H.-T., Eberl T., Ruffert M., Fryer C.L., 1999, ApJ 527, 

L39 

Kalogera V., Lorimer D.R., 2000, ApJ 530, 890 
Kochanek C.S., 1992, ApJ 398, 234 
Lai D., 1994, MNRAS 270, 611 

Shapiro S.L., 1995, ApJ 443, 705 



73 



Lai D. 
Lai D. 
Lai D. 
Lai D. 



Wiseman A.G., 1996, Phys. Rev. D 54, 3958 
Rasio F.A., Shapiro S.L., 1994a, ApJ 420, 811 
Rasio F.A., Shapiro S.L., 1994b, ApJ 423, 344 
Landau L.D., Lifschitz E.M., 1991, Lehrbuch der Theoretischen 
Physik, Band VI (Hydrodynamik, 5. Auflage), Akademie 
Verlag, Berlin 
Lattimer J.M., Schramm D.N., 1974, ApJ 192, L145 
Lattimer J.M., Schramm D.N., 1976, ApJ 210, 549 
Lattimer J.M., Swesty F.D., 1991, Nucl. Phys. A 535, 331 
Lattimer J.M., Mackie F., Ravenhall D.G., Schramm D.N., 

1977, ApJ 213, 225 
Lee H.K., Wijers R.A.M.J., Brown G.E., 2000, Physics 

Reports 325, 83 
Li L.-X., 2000, ApJ 533, L115 
Li L.-X., Paczyhski B., 2000, ApJ 534, L197 
Lombardi J.C., Rasio F.A., Shapiro S.L., 1997, Phys. 

Rev. D 56, 3416 
Marronetti P., Mathews G.J., Wilson J.R., 1999, 

Phys.Rev. D 60, 087301 
Meszaros P., Rees M.J., 1992, ApJ 397, 570 
Meszaros P., Rees M.J., 1997, ApJ 482, L29 
Meyer B.S., 1989, ApJ 343, 254 



1991, Prog. Theor. Phys. 86, 
, 1997, ApJ 490, 311 
1990, Prog. Theor. Phys. 83, 906 
1999, Prog. Theor. Phys. Suppl. 136, 



257 



Qian Y.-Z. 
Rasio F.A. 
Rasio F.A. 
Rasio F.A. 
Rosswog S., 



Ruffert M. 
Ruffert M. 
Ruffert M. 
Ruffert M 



ApJ 518 
1996, ApJ 471, 331 
1992, ApJ 401, 226 

1994, ApJ 432, 242 

1995, ApJ 438, 887 
Thielemann F.-K., Piran T 



356 



2000, 



Nakamura T., Oohara K., 
New K.C.B., Tohline J.E.^ 
Oohara K., Nakamura T., 
Oohara K., Nakamura T., 
270 

Paczyriski B., 1991, Acta Astronomica 41, 
Paczyhski B., Xu G., 1994, ApJ 427, 708 
Popham R., Woosley S.E., Fryer C, 1999, 
Woosley S.E., 
Shapiro S.L. 
Shapiro S.L. 
Shapiro S.L. 
Davies M.B. 
A&A 360, 171 

Rosswog S., Liebendorfer M., Thielemann F.-K., Davies M.B., 

Benz W., Piran T., 1999, A&A 341, 499 
Ruffert M., 1992, A&A 265, 82 

Janka H.-T., 1998, A&A 338, 535 
Janka H.-T., 1999, A&A 344, 573 
Janka H.-T., Schafer G., 1996, A&A 311, 532 
Janka H.-T., Takahashi K., Schafer G., 1997a, 
A&A 319, 122 

Ruffert M., Rampp M., Janka H.-T., 1997b, A&A 321, 991 
Sawyer R.F., 1980, ApJ 237, 187 

Shapiro S.L., Teukolsky S.A., 1983, Black Holes, White Dwarfs, 

and Neutron Stars, John Wiley & Sons, New York 
Shen H., Toki H., Oyamatsu K., Sumiyoshi K., 1998, Nucl. 

Phys. A 637, 435 
Shibata M., 1999, Phys. Rev. D 60, 104052 
Shibata M., Uryu K., 2000, Phys. Rev. D 61, 064001 
Shibata M., Uryri K., 2001, in Proc. 20th Texas Symposium on 

Relativistic Astrophysic s and Cosmology, cd. H. Martel & 

J.C. Wheeler, in press ( |istro-ph/010440c| ) 
Shibata M., Baumgarte T.W., Shapiro 

Phys.Rev. D 58, 023002 
Shibata M., Baumgarte T.W., Shapiro 

Phys.Rev. D 61, 044012 
Sumiyoshi K., Yamada S., Suzuki H., Hillebrandt W., 

A&A 334, 159 

Swesty F.D., Wang E.Y.M., Calder A.C., 2000, ApJ 541, 937 
Symbalisty E., Schramm D.N., 1982, Astrophys. Lett. 22, 143 
Taniguchi K., Nakamura T., 1996, Progress Theor. Phys. 96, 
693 

Thorne K.S., 1974, ApJ 191, 507 

Thorne K.S., 1995, in Procs. of the Snowmass 95 Summer 
Study on Particle and Nuclear Astrophysics and 
Cosmology, ed. E.W. Kolb & R. Peccei, World Scientific, 
Singapore, p. 398 

Thorne K.S., 1998, Phys. Rev. D 58, 124031 

Tscharnuter W.M., Winkler K.-H., 1979, Comp. Phys. 
Comm. 18, 171 

van den Horn L.J., van Weert C.G., 1981, Physics Lett. 84A, 
226 

Waxman E., BahcaU J.N., 1997, Phys. Rev. Lett. 78, 2292 

Waxman E., Bahcall J.N., 2000, ApJ 541, 707 

Weber F., 1999, Pulsars as Astrophysical Laboratories for 

Nuclear and Particle Physics, TOP, Bristol 
Woosley S.E., 1993a, ApJ 405, 273 
Woosley S.E., 1993b, A&AS 97, 205 

Zhuge X., Centrella J.M., McMillan S.L.W., 1994, Phys. 

Rev. D 50, 6247 
Zhuge X., Centrella J.M., McMillan S.L.W., 1996, Phys. 

Rev. D 54, 7261 



S.L. 



S.L. 



1998, 



2000, 



1998, 



This figure "setl.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "set2.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setS.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "set4.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setS.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "set6.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setT.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setS.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "set9.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setlO.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setll.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setl2.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



This figure "setlS.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/0106229v2 



